Gam models: summary

Birthweight

2007-2020

m_BW_bam5_2 <- bam(BW ~ s(GA_weeks) +  s(mat_age, k=7) + s(birth_Y_M_num) + s(month) + s(mean_ssep2) + parity_cat + sex + Urban3 + Language + mother_nationality_cat2, data=bevn_eco_in7)  
  summary(m_BW_bam5_2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## BW ~ s(GA_weeks) + s(mat_age, k = 7) + s(birth_Y_M_num) + s(month) + 
##     s(mean_ssep2) + parity_cat + sex + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                       Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)                                          3335.2704     0.8361 3989.044  < 2e-16 ***
## parity_cat(1,2]                                       127.4692     0.8344  152.762  < 2e-16 ***
## parity_cat(2,3]                                       173.6260     1.2772  135.939  < 2e-16 ***
## parity_cat(3,20]                                      199.6977     2.2149   90.160  < 2e-16 ***
## sex2                                                 -147.0101     0.7458 -197.128  < 2e-16 ***
## Urban3                                                  0.6608     0.7906    0.836    0.403    
## LanguageFrench                                        -48.0924     0.9182  -52.374  < 2e-16 ***
## LanguageItalian                                       -63.1878     2.0493  -30.834  < 2e-16 ***
## mother_nationality_cat2Africa                          11.0168     2.2735    4.846 1.26e-06 ***
## mother_nationality_cat2Asia                           -17.9253     2.0596   -8.703  < 2e-16 ***
## mother_nationality_cat2Europe                          32.5679     0.8501   38.310  < 2e-16 ***
## mother_nationality_cat2Northern America                65.0093     4.8546   13.391  < 2e-16 ***
## mother_nationality_cat2Oceania                         63.4678    11.1288    5.703 1.18e-08 ***
## mother_nationality_cat2Southern and Central America    47.0493     2.7991   16.809  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df         F p-value    
## s(GA_weeks)      8.826  8.988 73217.048  <2e-16 ***
## s(mat_age)       5.178  5.693    20.286  <2e-16 ***
## s(birth_Y_M_num) 8.285  8.859    39.009  <2e-16 ***
## s(month)         5.893  7.053     3.183  0.0022 ** 
## s(mean_ssep2)    7.634  8.534    43.693  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.401   Deviance explained = 40.1%
## fREML = 8.1205e+06  Scale est. = 1.5259e+05  n = 1099328
  gam.check(m_BW_bam5_2) 

## 
## Method: fREML   Optimizer: perf newton
## full convergence after 11 iterations.
## Gradient range [-0.07345475,0.05455442]
## (score 8120451 & scale 152594.7).
## Hessian positive definite, eigenvalue range [1.608974,549654.6].
## Model rank =  56 / 56 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                    k'  edf k-index p-value  
## s(GA_weeks)      9.00 8.83    0.98   0.095 .
## s(mat_age)       6.00 5.18    0.99   0.325  
## s(birth_Y_M_num) 9.00 8.28    0.98   0.095 .
## s(month)         9.00 5.89    1.02   0.965  
## s(mean_ssep2)    9.00 7.63    1.03   0.995  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  plot(m_BW_bam5_2, select=1, ylim=c(0, 4000),shade = TRUE, shade.col = "lightblue",  shift = coef(m_BW_bam5_2)[1])

  plot(m_BW_bam5_2, select=2, ylim=c(3200, 3400),shade = TRUE, shade.col = "lightblue",  shift = coef(m_BW_bam5_2)[1])

  plot(m_BW_bam5_2, select=3, ylim=c(3300, 3380),shade = TRUE, shade.col = "lightblue",  shift = coef(m_BW_bam5_2)[1]  ,  xlab  = "1 : 2007-01,     25: 2009-01,     50: 2011-02,    75: 2013-03,   100: 2015-04,    125: 2017-05, 150: 2019-06")

  plot(m_BW_bam5_2, select=4, ylim=c(3325, 3345),shade = TRUE, shade.col = "lightblue",  shift = coef(m_BW_bam5_2)[1])

  plot(m_BW_bam5_2, select=5, ylim=c(3200, 3400),shade = TRUE, shade.col = "lightblue",  shift = coef(m_BW_bam5_2)[1])

Adjusting for the effect of the huge dataset (n=1,099,368)

Cohen’s d formula:

\(d=\frac{beta}{SD}\) and \(SD=\sqrt{n}*se\)

\(d=\frac{beta}{\sqrt{n}*se}\)

  • beta: estimate
  • se: standard error
  • n: population/observation size
  • SD: standard deviation
#getting the CI 95% and d (interpret only for the non-smooth parameters)
  beta <- coef(m_BW_bam5_2)
  Vb <- vcov(m_BW_bam5_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  n <- 1099328
  d <- (abs(beta)/(sqrt(n)*se))
  my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))
  round(head(my.ci,14), 2)
##                                                         beta     lci     uci    d
## (Intercept)                                          3335.27 3333.63 3336.91 3.80
## parity_cat(1,2]                                       127.47  125.83  129.10 0.15
## parity_cat(2,3]                                       173.63  171.12  176.13 0.13
## parity_cat(3,20]                                      199.70  195.36  204.04 0.09
## sex2                                                 -147.01 -148.47 -145.55 0.19
## Urban3                                                  0.66   -0.89    2.21 0.00
## LanguageFrench                                        -48.09  -49.89  -46.29 0.05
## LanguageItalian                                       -63.19  -67.20  -59.17 0.03
## mother_nationality_cat2Africa                          11.02    6.56   15.47 0.00
## mother_nationality_cat2Asia                           -17.93  -21.96  -13.89 0.01
## mother_nationality_cat2Europe                          32.57   30.90   34.23 0.04
## mother_nationality_cat2Northern America                65.01   55.49   74.52 0.01
## mother_nationality_cat2Oceania                         63.47   41.66   85.28 0.01
## mother_nationality_cat2Southern and Central America    47.05   41.56   52.54 0.02

Interpretation of m_BW_bam5_2 model :

Why choosing m_BW_bam5_2? It has the lowest AIC and also lower k value for maternal age, which seems more appropriate by looking at the graph.

  • Non-smooth terms

Increasing parity increases birthweight compared to first births: from (126g, 129g) at second parity, (171g, 176g) at third parity, to (195g, 204g) for parities over 3.

Being born female decreases by (-148g, -146g) compared to males.

Urbanicity has no effect on birthweight.

Language: Compared to babies born in German(+Romansh)-speaking Switzerland, being born in the French-speaking part decreases BW by (-50g,-46g), while being born in the Italian-speaking part decreases BW by (-67g,-59g).

Maternal nationality: compared to babies born from Swiss mothers, babies born from Asian mothers are lighter (-22g, -14g). Babies born from other-nationalities mothers are heavier: (7g, 16g) for Africa, (31g, 34g) for other European countries, (55g, 75g) for Northern America, (42g, 86g) for Oceania, and (42g-53g) for Southern/Central America.

  • Smooth terms

GA: Obviously, birthweight increases with increasing gestational age, with a plateau at around 800g for babies between 20 and 25 gestational weeks, and another plateau at around 4000g for babies above 45 weeks.

Maternal age: birthweight increases non-linearly with maternal age: for younger mothers (<25yo) it increases almost linearly from 3280g to 3350g, then declines very slightly until around 3340g when the mother is 25-40yo. Older women have heavier babies, increasing from 3340g to around 3360g for women 40-50yo (keep in mind we have excluded a few mothers above 50 - improbable).

Time (month+year combined): globally, birthweight is decreasing between 2007 and 2019 from 3350 to 3330 which is a minor change, and then increases from 2019 to reach 3350g. There are some other minor variations throughtout the whole time but this doesnt seem very relevant: only 20g variation.

Seasonal/month: minor sinusoidal seasonal variation between 3330-3340g, with lower birthweight in March and August. These difference dont seem much relevant either (10g variation will probably not impact the baby’s health).

Mean ssep: Birthweight increases by increasing ssep, linearly between ssep=25-45 approx (3250-3340g), and then there is some minor birthweight increase for higher ssep, until reaching a BW of around 3360g. Overall, a +100g increase is probably worht mentioning.

2008-2010

Because birth_Y_M_num and month variables correlated too much, moth variable was not kept in the model. Mean ssep was set as linear variables (when included as a smooth, the graph clearly showed a linear association).

  m_BW_bam_2009_2 <- bam(BW ~ s(GA_weeks) +  s(mat_age, k=7) + s(birth_Y_M_num) + mean_ssep2 + parity_cat + sex + Urban3  + Language + mother_nationality_cat2, data=bevn_eco_in7_2008_10)
  summary(m_BW_bam_2009_2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## BW ~ s(GA_weeks) + s(mat_age, k = 7) + s(birth_Y_M_num) + mean_ssep2 + 
##     parity_cat + sex + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                                          3278.4522     6.4015 512.135  < 2e-16 ***
## mean_ssep2                                              0.8785     0.1043   8.419  < 2e-16 ***
## parity_cat(1,2]                                       122.1534     1.8847  64.813  < 2e-16 ***
## parity_cat(2,3]                                       172.2947     2.8838  59.746  < 2e-16 ***
## parity_cat(3,20]                                      196.5618     4.9231  39.926  < 2e-16 ***
## sex2                                                 -146.2210     1.6800 -87.037  < 2e-16 ***
## Urban3                                                  0.8197     1.7447   0.470  0.63846    
## LanguageFrench                                        -45.9768     2.0608 -22.310  < 2e-16 ***
## LanguageItalian                                       -65.0958     4.3465 -14.977  < 2e-16 ***
## mother_nationality_cat2Africa                          41.2721     6.0377   6.836 8.18e-12 ***
## mother_nationality_cat2Asia                           -15.2504     4.9412  -3.086  0.00203 ** 
## mother_nationality_cat2Europe                          36.4670     1.9396  18.801  < 2e-16 ***
## mother_nationality_cat2Northern America                82.8570    10.6237   7.799 6.25e-15 ***
## mother_nationality_cat2Oceania                         54.3020    24.9696   2.175  0.02965 *  
## mother_nationality_cat2Southern and Central America    60.6874     6.0057  10.105  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df         F  p-value    
## s(GA_weeks)      8.246  8.808 14713.688  < 2e-16 ***
## s(mat_age)       3.134  3.771     7.669 1.06e-05 ***
## s(birth_Y_M_num) 2.838  3.533    14.969  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.397   Deviance explained = 39.7%
## fREML = 1.6251e+06  Scale est. = 1.548e+05  n = 219790
   plot(m_BW_bam_2009_2, select=1, ylim=c(0, 4000),shade = TRUE, shade.col = "lightblue",shift = coef(m_BW_bam_2009_2)[1]) #GA

  plot(m_BW_bam_2009_2, select=2, ylim=c(3200, 3350),shade = TRUE, shade.col = "lightblue",shift = coef(m_BW_bam_2009_2)[1]) #mat age

  plot(m_BW_bam_2009_2, select=3, ylim=c(3250, 3350),shade = TRUE, shade.col = "lightblue", shift = coef(m_BW_bam_2009_2)[1],  xlab  = "13 : 2008-01,   20: 2008-08,   30: 2009-06,    40: 2010-04,   50: 2011-02") # time

  gam.check(m_BW_bam_2009_2)

## 
## Method: fREML   Optimizer: perf newton
## full convergence after 9 iterations.
## Gradient range [-0.002922242,0.001533332]
## (score 1625083 & scale 154802.2).
## Hessian positive definite, eigenvalue range [0.2911485,109886].
## Model rank =  39 / 39 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                    k'  edf k-index p-value  
## s(GA_weeks)      9.00 8.25    0.99    0.38  
## s(mat_age)       6.00 3.13    0.98    0.08 .
## s(birth_Y_M_num) 9.00 2.84    0.99    0.36  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  concurvity(m_BW_bam_2009_2, full=TRUE) #values are now fine.
##               para s(GA_weeks) s(mat_age) s(birth_Y_M_num)
## worst    0.9828236 0.015759637 0.13707246      0.005557451
## observed 0.9828236 0.006364574 0.02758171      0.005155801
## estimate 0.9828236 0.009119694 0.08566593      0.003170851
  # Summary of estimates, CI95% and Cohen's d
  beta <- coef(m_BW_bam_2009_2)
  Vb <- vcov(m_BW_bam_2009_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  n <- 219790
  d <- (abs(beta)/(sqrt(n)*se))
  my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))
  round(head(my.ci,15), 2)  
##                                                         beta     lci     uci    d
## (Intercept)                                          3278.45 3265.90 3291.00 1.09
## mean_ssep2                                              0.88    0.67    1.08 0.02
## parity_cat(1,2]                                       122.15  118.46  125.85 0.14
## parity_cat(2,3]                                       172.29  166.64  177.95 0.13
## parity_cat(3,20]                                      196.56  186.91  206.21 0.09
## sex2                                                 -146.22 -149.51 -142.93 0.19
## Urban3                                                  0.82   -2.60    4.24 0.00
## LanguageFrench                                        -45.98  -50.02  -41.94 0.05
## LanguageItalian                                       -65.10  -73.61  -56.58 0.03
## mother_nationality_cat2Africa                          41.27   29.44   53.11 0.01
## mother_nationality_cat2Asia                           -15.25  -24.94   -5.57 0.01
## mother_nationality_cat2Europe                          36.47   32.66   40.27 0.04
## mother_nationality_cat2Northern America                82.86   62.03  103.68 0.02
## mother_nationality_cat2Oceania                         54.30    5.36  103.24 0.00
## mother_nationality_cat2Southern and Central America    60.69   48.92   72.46 0.02

2017-2020

We also did not use the month variable because of correlation between month and time (birth_Y_M_num) variables.

  m_BW_bam_covid_2 <- bam(BW ~ s(GA_weeks) +  s(mat_age, k=7) + s(birth_Y_M_num) + s(mean_ssep2) + parity_cat + sex + Urban3  + Language + mother_nationality_cat2, data=bevn_eco_in7_2017_20)
  summary(m_BW_bam_covid_2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## BW ~ s(GA_weeks) + s(mat_age, k = 7) + s(birth_Y_M_num) + s(mean_ssep2) + 
##     parity_cat + sex + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                        Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)                                          3341.78975    1.51250 2209.446  < 2e-16 ***
## parity_cat(1,2]                                       132.59338    1.50324   88.205  < 2e-16 ***
## parity_cat(2,3]                                       180.82725    2.29301   78.860  < 2e-16 ***
## parity_cat(3,20]                                      202.25754    4.00225   50.536  < 2e-16 ***
## sex2                                                 -144.71388    1.34661 -107.465  < 2e-16 ***
## Urban3                                                 -0.09571    1.42256   -0.067  0.94636    
## LanguageFrench                                        -51.18569    1.63285  -31.348  < 2e-16 ***
## LanguageItalian                                       -62.48479    3.94863  -15.824  < 2e-16 ***
## mother_nationality_cat2Africa                          -6.51636    3.76791   -1.729  0.08373 .  
## mother_nationality_cat2Asia                           -16.40941    3.54593   -4.628 3.70e-06 ***
## mother_nationality_cat2Europe                          27.23301    1.52212   17.892  < 2e-16 ***
## mother_nationality_cat2Northern America                59.93225    9.02825    6.638 3.18e-11 ***
## mother_nationality_cat2Oceania                         60.90901   22.28696    2.733  0.00628 ** 
## mother_nationality_cat2Southern and Central America    38.21351    5.36683    7.120 1.08e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df         F p-value    
## s(GA_weeks)      8.504  8.906 22216.026  <2e-16 ***
## s(mat_age)       4.935  5.500    23.982  <2e-16 ***
## s(birth_Y_M_num) 5.636  6.787     9.851  <2e-16 ***
## s(mean_ssep2)    5.067  6.235    15.020  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.399   Deviance explained = 39.9%
## fREML = 2.4557e+06  Scale est. = 1.5057e+05  n = 332752
  plot(m_BW_bam_covid_2, select=1, ylim=c(0, 4000),shade = TRUE, shade.col = "lightblue",shift = coef(m_BW_bam_covid_2)[1])

  plot(m_BW_bam_covid_2, select=2, ylim=c(3200, 3400),shade = TRUE, shade.col = "lightblue",shift = coef(m_BW_bam_covid_2)[1])

  plot(m_BW_bam_covid_2, select=3, ylim=c(3310, 3360),shade = TRUE, shade.col = "lightblue", shift = coef(m_BW_bam_covid_2)[1],xlab  = "120 : 2016-12,   130: 2017-10,   140: 2018-08,   150: 2019-06,   160: 2020-04")

  plot(m_BW_bam_covid_2, select=4, ylim=c(3200, 3400),shade = TRUE, shade.col = "lightblue", shift = coef(m_BW_bam_covid_2)[1])

  gam.check(m_BW_bam_covid_2)

## 
## Method: fREML   Optimizer: perf newton
## full convergence after 12 iterations.
## Gradient range [-0.02069349,0.01034373]
## (score 2455716 & scale 150568.5).
## Hessian positive definite, eigenvalue range [0.9600157,166367].
## Model rank =  47 / 47 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                    k'  edf k-index p-value
## s(GA_weeks)      9.00 8.50    1.00    0.39
## s(mat_age)       6.00 4.93    1.03    0.99
## s(birth_Y_M_num) 9.00 5.64    1.00    0.39
## s(mean_ssep2)    9.00 5.07    1.00    0.54
  concurvity(m_BW_bam_covid_2, full=TRUE) #values are now fine.
##               para s(GA_weeks) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst    0.8024701 0.018706682 0.10865948     0.0012092151    0.17776696
## observed 0.8024701 0.007194495 0.06217053     0.0006631716    0.08814200
## estimate 0.8024701 0.010875530 0.06838436     0.0009362815    0.08021785
  # Summary of estimates, CI95% and Cohen's d
  beta <- coef(m_BW_bam_covid_2)
  Vb <- vcov(m_BW_bam_covid_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  n <- 332752
  my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))
## Warning in cbind(beta, lci = beta - 1.96 * se, uci = beta + 1.96 * se, d): number of rows of result
## is not a multiple of vector length (arg 4)
  round(head(my.ci,14), 2)
##                                                         beta     lci     uci    d
## (Intercept)                                          3341.79 3338.83 3344.75 1.09
## parity_cat(1,2]                                       132.59  129.65  135.54 0.02
## parity_cat(2,3]                                       180.83  176.33  185.32 0.14
## parity_cat(3,20]                                      202.26  194.41  210.10 0.13
## sex2                                                 -144.71 -147.35 -142.07 0.09
## Urban3                                                 -0.10   -2.89    2.70 0.19
## LanguageFrench                                        -51.19  -54.39  -47.98 0.00
## LanguageItalian                                       -62.48  -70.23  -54.74 0.05
## mother_nationality_cat2Africa                          -6.52  -13.90    0.87 0.03
## mother_nationality_cat2Asia                           -16.41  -23.36   -9.46 0.01
## mother_nationality_cat2Europe                          27.23   24.25   30.22 0.01
## mother_nationality_cat2Northern America                59.93   42.24   77.63 0.04
## mother_nationality_cat2Oceania                         60.91   17.23  104.59 0.02
## mother_nationality_cat2Southern and Central America    38.21   27.69   48.73 0.00

Stillbirth

2007-2020

m_SB_bam2 <- bam(stillbirth~s(GA_weeks)+ s(mat_age, k=7) + s(month) + birth_Y_M_num + mean_ssep2 + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in6)
  summary(m_SB_bam2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## stillbirth ~ s(GA_weeks) + s(mat_age, k = 7) + s(month) + birth_Y_M_num + 
##     mean_ssep2 + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                          -6.619e+00  1.456e-01 -45.452  < 2e-16 ***
## birth_Y_M_num                                         4.545e-05  3.972e-04   0.114  0.90891    
## mean_ssep2                                           -6.094e-03  2.332e-03  -2.614  0.00896 ** 
## Urban3                                                3.814e-02  3.892e-02   0.980  0.32712    
## LanguageFrench                                        8.110e-02  4.404e-02   1.842  0.06554 .  
## LanguageItalian                                      -2.043e-02  1.024e-01  -0.200  0.84185    
## mother_nationality_cat2Africa                         2.275e-01  9.350e-02   2.434  0.01495 *  
## mother_nationality_cat2Asia                          -3.683e-02  1.016e-01  -0.362  0.71708    
## mother_nationality_cat2Europe                         3.175e-02  4.253e-02   0.747  0.45532    
## mother_nationality_cat2Northern America              -3.488e-01  2.851e-01  -1.223  0.22121    
## mother_nationality_cat2Oceania                       -1.733e-01  7.207e-01  -0.241  0.80993    
## mother_nationality_cat2Southern and Central America  -9.889e-02  1.351e-01  -0.732  0.46426    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df    Chi.sq p-value    
## s(GA_weeks) 7.589  8.143 16734.305  <2e-16 ***
## s(mat_age)  3.147  3.820    13.061  0.0105 *  
## s(month)    2.627  3.267     7.216  0.0701 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.133   Deviance explained = 32.1%
## fREML = 1.5655e+06  Scale est. = 1         n = 1102617
  plot(m_SB_bam2, trans = plogis, select=1, shift = coef(m_SB_bam2)[1], shade = TRUE, shade.col = "lightblue") #GA

  plot(m_SB_bam2, trans = plogis, select=2, ylim=c(0, 0.003), shift = coef(m_SB_bam2)[1], shade = TRUE, shade.col = "lightblue") #mat age

  plot(m_SB_bam2, trans = plogis, select=3, ylim=c(0.0005, 0.0015), shift = coef(m_SB_bam2)[1], shade = TRUE, shade.col = "lightblue") #month

  k.check(m_SB_bam2) 
##             k'      edf   k-index p-value
## s(GA_weeks)  9 7.588646 0.9509095  0.6175
## s(mat_age)   6 3.146678 0.9511237  0.6800
## s(month)     9 2.626666 0.9507136  0.5325
  concurvity(m_SB_bam2)
##               para s(GA_weeks) s(mat_age)     s(month)
## worst    0.9829412 0.010976422 0.05328595 0.0045643900
## observed 0.9829412 0.005611189 0.01081304 0.0005238286
## estimate 0.9829412 0.007092571 0.03420192 0.0034099344

Calculation of Cohen’s d for OR:

\(d= \log(OR) * \frac{\sqrt{3}}{\pi}\)

\(d = \frac{log(OR)}{1.81}\)

#table to show OR, confidence interval and Cohen's d
  OR <- exp(coef(m_SB_bam2))
  beta <- coef(m_SB_bam2)
  Vb <- vcov(m_SB_bam2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 12), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.00 0.00 0.00 3.66
## birth_Y_M_num                                        1.00 1.00 1.00 0.00
## mean_ssep2                                           0.99 0.99 1.00 0.00
## Urban3                                               1.04 0.96 1.12 0.02
## LanguageFrench                                       1.08 0.99 1.18 0.04
## LanguageItalian                                      0.98 0.80 1.20 0.01
## mother_nationality_cat2Africa                        1.26 1.05 1.51 0.13
## mother_nationality_cat2Asia                          0.96 0.79 1.18 0.02
## mother_nationality_cat2Europe                        1.03 0.95 1.12 0.02
## mother_nationality_cat2Northern America              0.71 0.40 1.23 0.19
## mother_nationality_cat2Oceania                       0.84 0.20 3.45 0.10
## mother_nationality_cat2Southern and Central America  0.91 0.70 1.18 0.05

Time variables (birth_Y_M_num) and mean ssep variable were put as linear terms, because the model including them as smoothed terms showed that the relationship between these and stillbirth was linear.

#Adding BW –> discussing this.

 m_SB_bam3 = update(m_SB_bam2, . ~ . + s(BW))
plot(m_SB_bam3, trans = plogis, select=4, ylim=c(0, 1), shift = coef(m_SB_bam3)[1], shade = TRUE, shade.col = "lightblue") #BW with whole CI

  plot(m_SB_bam3, trans = plogis, select=4, ylim=c(0, 0.15), shift = coef(m_SB_bam3)[1], shade = TRUE, shade.col = "lightblue") #BW

  concurvity(m_SB_bam3)
##               para s(GA_weeks) s(mat_age)     s(month)     s(BW)
## worst    0.9829455   0.8787415 0.05483844 0.0045794820 0.8787530
## observed 0.9829455   0.4617487 0.01746103 0.0006330721 0.5160650
## estimate 0.9829455   0.4140313 0.03533273 0.0034246658 0.3330819
  k.check(m_SB_bam3)
##             k'      edf   k-index p-value
## s(GA_weeks)  9 7.293198 0.9949575  1.0000
## s(mat_age)   6 2.760090 0.9692146  0.5475
## s(month)     9 2.699310 0.9677931  0.2725
## s(BW)        9 6.933464 0.9951825  1.0000

When we add birthweight, there is too much correlation betwwen gestational age and birthweight (concurvity shows 0.88 which is to close to 1): it is not a good idea to include BW in the model. But we can have a look at the smoothed curve for stillbirth depending on birthweight: it has a clear U shaped, with high stillbirth risk (2-11%) for very low birthweight babies (<2kg), and moderately high stillbirth risk (5-8%) for very high birthweight babies (>8kg).

2008-2010

m_SB_bam_2009_2<-bam(stillbirth~s(GA_weeks)+ s(mat_age, k=7) + birth_Y_M_num + mean_ssep2 + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in6_2008_10)
summary(m_SB_bam_2009_2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## stillbirth ~ s(GA_weeks) + s(mat_age, k = 7) + birth_Y_M_num + 
##     mean_ssep2 + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                          -6.408576   0.339228 -18.892   <2e-16 ***
## birth_Y_M_num                                        -0.003879   0.004052  -0.957   0.3385    
## mean_ssep2                                           -0.007007   0.005230  -1.340   0.1803    
## Urban3                                                0.174863   0.086722   2.016   0.0438 *  
## LanguageFrench                                        0.033909   0.100425   0.338   0.7356    
## LanguageItalian                                       0.043843   0.208702   0.210   0.8336    
## mother_nationality_cat2Africa                        -0.307892   0.262929  -1.171   0.2416    
## mother_nationality_cat2Asia                          -0.623295   0.292794  -2.129   0.0333 *  
## mother_nationality_cat2Europe                         0.065891   0.093736   0.703   0.4821    
## mother_nationality_cat2Northern America              -0.832096   0.740930  -1.123   0.2614    
## mother_nationality_cat2Oceania                       -8.197695  52.128818  -0.157   0.8750    
## mother_nationality_cat2Southern and Central America  -0.318662   0.314474  -1.013   0.3109    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df  Chi.sq p-value    
## s(GA_weeks) 5.937  6.945 3434.41  <2e-16 ***
## s(mat_age)  2.034  2.551    3.07   0.301    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.149   Deviance explained = 32.5%
## fREML = 3.1243e+05  Scale est. = 1         n = 220464
   plot(m_SB_bam_2009_2, trans = plogis, select=1, ylim=c(0, 0.9) ,shift = coef(m_SB_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #GA

  plot(m_SB_bam_2009_2, trans = plogis, select=2, ylim=c(0, 0.005), shift = coef(m_SB_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #mat age

  k.check(m_SB_bam_2009_2)
##             k'      edf   k-index p-value
## s(GA_weeks)  9 5.936556 0.8528645  0.0150
## s(mat_age)   6 2.033857 0.9027108  0.2975
  concurvity(m_SB_bam_2009_2)
##               para s(GA_weeks)  s(mat_age)
## worst    0.9840551 0.010398851 0.061160016
## observed 0.9840551 0.005436465 0.005779778
## estimate 0.9840551 0.006773190 0.035691703
  OR <- exp(coef(m_SB_bam_2009_2))
  beta <- coef(m_SB_bam_2009_2)
  Vb <- vcov(m_SB_bam_2009_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 12), 2)
##                                                        OR  lci          uci    d
## (Intercept)                                          0.00 0.00 0.000000e+00 3.54
## birth_Y_M_num                                        1.00 0.99 1.000000e+00 0.00
## mean_ssep2                                           0.99 0.98 1.000000e+00 0.00
## Urban3                                               1.19 1.00 1.410000e+00 0.10
## LanguageFrench                                       1.03 0.85 1.260000e+00 0.02
## LanguageItalian                                      1.04 0.69 1.570000e+00 0.02
## mother_nationality_cat2Africa                        0.73 0.44 1.230000e+00 0.17
## mother_nationality_cat2Asia                          0.54 0.30 9.500000e-01 0.34
## mother_nationality_cat2Europe                        1.07 0.89 1.280000e+00 0.04
## mother_nationality_cat2Northern America              0.44 0.10 1.860000e+00 0.46
## mother_nationality_cat2Oceania                       0.00 0.00 6.497288e+40 4.53
## mother_nationality_cat2Southern and Central America  0.73 0.39 1.350000e+00 0.18

As for the previous model (including all timing 2007-2020), we kept time(birth_Y_M_num) and mean_ssep2 as linear terms.

2017-2020

m_SB_bam_covid_2 <-bam(stillbirth~s(GA_weeks)+ mat_age + birth_Y_M_num +  mean_ssep2 + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in6_2017_20)
summary(m_SB_bam_covid_2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## stillbirth ~ s(GA_weeks) + mat_age + birth_Y_M_num + mean_ssep2 + 
##     Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                          -6.794848   0.471639 -14.407  < 2e-16 ***
## mat_age                                               0.007085   0.006871   1.031 0.302447    
## birth_Y_M_num                                        -0.002495   0.002496  -1.000 0.317501    
## mean_ssep2                                           -0.001567   0.004273  -0.367 0.713770    
## Urban3                                               -0.121121   0.071465  -1.695 0.090107 .  
## LanguageFrench                                        0.168834   0.078861   2.141 0.032282 *  
## LanguageItalian                                      -0.244986   0.230082  -1.065 0.286978    
## mother_nationality_cat2Africa                         0.552349   0.149243   3.701 0.000215 ***
## mother_nationality_cat2Asia                           0.192085   0.167090   1.150 0.250312    
## mother_nationality_cat2Europe                         0.022826   0.078368   0.291 0.770843    
## mother_nationality_cat2Northern America               0.274686   0.439897   0.624 0.532344    
## mother_nationality_cat2Oceania                       -8.133629  46.445173  -0.175 0.860983    
## mother_nationality_cat2Southern and Central America   0.044082   0.242353   0.182 0.855667    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq p-value    
## s(GA_weeks) 7.421  8.111   4940  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.127   Deviance explained = 33.1%
## fREML = 4.7434e+05  Scale est. = 1         n = 333721
 plot(m_SB_bam_covid_2, trans = plogis, select=1, ylim=c(0, 1) ,shift = coef(m_SB_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue") #GA

  OR <- exp(coef(m_SB_bam_covid_2))
  beta <- coef(m_SB_bam_covid_2)
  Vb <- vcov(m_SB_bam_covid_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 13), 2)
##                                                        OR  lci          uci    d
## (Intercept)                                          0.00 0.00 0.000000e+00 3.75
## mat_age                                              1.01 0.99 1.020000e+00 0.00
## birth_Y_M_num                                        1.00 0.99 1.000000e+00 0.00
## mean_ssep2                                           1.00 0.99 1.010000e+00 0.00
## Urban3                                               0.89 0.77 1.020000e+00 0.07
## LanguageFrench                                       1.18 1.01 1.380000e+00 0.09
## LanguageItalian                                      0.78 0.50 1.230000e+00 0.14
## mother_nationality_cat2Africa                        1.74 1.30 2.330000e+00 0.31
## mother_nationality_cat2Asia                          1.21 0.87 1.680000e+00 0.11
## mother_nationality_cat2Europe                        1.02 0.88 1.190000e+00 0.01
## mother_nationality_cat2Northern America              1.32 0.56 3.120000e+00 0.15
## mother_nationality_cat2Oceania                       0.00 0.00 1.005865e+36 4.49
## mother_nationality_cat2Southern and Central America  1.05 0.65 1.680000e+00 0.02
 k.check(m_SB_bam_covid_2)
##             k'      edf   k-index p-value
## s(GA_weeks)  9 7.420665 0.8755089    0.02
 concurvity(m_SB_bam_2009_2)
##               para s(GA_weeks)  s(mat_age)
## worst    0.9840551 0.010398851 0.061160016
## observed 0.9840551 0.005436465 0.005779778
## estimate 0.9840551 0.006773190 0.035691703

Why we chose this model: Maternal age had edf values very close to one: we put it as a linear term instead of smoothed variable. As for the previous model (including all timing 2007-2020), we kept time (birth_Y_M_num) and mean_ssep2 variables as linear terms.

Gestational age in days (continuous)

2007-2020

m_GA_bam_days <- bam(GA_days ~ s(BW) +  s(mat_age) + s(birth_Y_M_num) + s(month) + s(mean_ssep2) + parity_cat + sex + Urban3 + mother_nationality_cat2 + Language, data=bevn_eco_in7)
summary(m_GA_bam_days)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GA_days ~ s(BW) + s(mat_age) + s(birth_Y_M_num) + s(month) + 
##     s(mean_ssep2) + parity_cat + sex + Urban3 + mother_nationality_cat2 + 
##     Language
## 
## Parametric coefficients:
##                                                       Estimate Std. Error   t value Pr(>|t|)    
## (Intercept)                                          276.04111    0.01752 15751.707  < 2e-16 ***
## parity_cat(1,2]                                       -1.99202    0.01754  -113.596  < 2e-16 ***
## parity_cat(2,3]                                       -2.49376    0.02687   -92.814  < 2e-16 ***
## parity_cat(3,20]                                      -2.56052    0.04653   -55.033  < 2e-16 ***
## sex2                                                   1.81258    0.01581   114.644  < 2e-16 ***
## Urban3                                                 0.09180    0.01657     5.539 3.04e-08 ***
## mother_nationality_cat2Africa                          0.60971    0.04764    12.798  < 2e-16 ***
## mother_nationality_cat2Asia                           -0.79664    0.04316   -18.460  < 2e-16 ***
## mother_nationality_cat2Europe                         -0.35922    0.01783   -20.147  < 2e-16 ***
## mother_nationality_cat2Northern America               -0.28620    0.10175    -2.813  0.00491 ** 
## mother_nationality_cat2Oceania                        -0.32052    0.23325    -1.374  0.16939    
## mother_nationality_cat2Southern and Central America   -1.59307    0.05865   -27.161  < 2e-16 ***
## LanguageFrench                                         0.10488    0.01928     5.441 5.29e-08 ***
## LanguageItalian                                       -0.20128    0.04297    -4.684 2.81e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df         F p-value    
## s(BW)            8.977  9.000 127518.85  <2e-16 ***
## s(mat_age)       8.044  8.641    347.60  <2e-16 ***
## s(birth_Y_M_num) 8.689  8.972    115.30  <2e-16 ***
## s(month)         6.958  8.043     12.76  <2e-16 ***
## s(mean_ssep2)    7.867  8.668     59.48  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.515   Deviance explained = 51.5%
## fREML = 3.8714e+06  Scale est. = 67.031    n = 1099328
  #plot(m_GA_bam_days, all.terms = TRUE,pages=1, shade = TRUE, shade.col = "lightblue")
  plot(m_GA_bam_days, select=1, shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_days)[1]) #BW

  plot(m_GA_bam_days, select=2, ylim=c(265, 280),shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_days)[1]) #mat age

  plot(m_GA_bam_days, select=3, ylim=c(274, 278),shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_days)[1],  xlab  = "1 : 2007-01,    25: 2009-01,     50: 2011-02,   75: 2013-03,     100: 2015-04,    125: 2017-05,   150: 2019-06") #birth_Y_M_num

  plot(m_GA_bam_days, select=4, ylim=c(275, 277),shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_days)[1]) #month

  plot(m_GA_bam_days, select=5, ylim=c(272, 280),shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_days)[1]) #ssep

  gam.check(m_GA_bam_days)

## 
## Method: fREML   Optimizer: perf newton
## full convergence after 13 iterations.
## Gradient range [-0.001200072,0.001194492]
## (score 3871436 & scale 67.03148).
## Hessian positive definite, eigenvalue range [1.310202,549654.5].
## Model rank =  59 / 59 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                    k'  edf k-index p-value  
## s(BW)            9.00 8.98    1.00   0.500  
## s(mat_age)       9.00 8.04    0.99   0.360  
## s(birth_Y_M_num) 9.00 8.69    0.98   0.085 .
## s(month)         9.00 6.96    1.01   0.815  
## s(mean_ssep2)    9.00 7.87    1.01   0.775  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  concurvity(m_GA_bam_days, full=TRUE)
##               para      s(BW) s(mat_age) s(birth_Y_M_num)    s(month) s(mean_ssep2)
## worst    0.8015268 0.04840718 0.12110433       0.03003088 0.028058758    0.17917903
## observed 0.8015268 0.01911969 0.07090416       0.01041702 0.001316444    0.04934839
## estimate 0.8015268 0.03133441 0.08072505       0.01100951 0.021033963    0.08954157
  #estimate, 95%CI and d
  beta <- coef(m_GA_bam_days)
  Vb <- vcov(m_GA_bam_days, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  n <- 1099328
  d <- (abs(beta)/(sqrt(n)*se))
  my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))
  round(head(my.ci,14), 2)
##                                                        beta    lci    uci     d
## (Intercept)                                          276.04 276.01 276.08 15.02
## parity_cat(1,2]                                       -1.99  -2.03  -1.96  0.11
## parity_cat(2,3]                                       -2.49  -2.55  -2.44  0.09
## parity_cat(3,20]                                      -2.56  -2.65  -2.47  0.05
## sex2                                                   1.81   1.78   1.84  0.11
## Urban3                                                 0.09   0.06   0.12  0.01
## mother_nationality_cat2Africa                          0.61   0.52   0.70  0.01
## mother_nationality_cat2Asia                           -0.80  -0.88  -0.71  0.02
## mother_nationality_cat2Europe                         -0.36  -0.39  -0.32  0.02
## mother_nationality_cat2Northern America               -0.29  -0.49  -0.09  0.00
## mother_nationality_cat2Oceania                        -0.32  -0.78   0.14  0.00
## mother_nationality_cat2Southern and Central America   -1.59  -1.71  -1.48  0.03
## LanguageFrench                                         0.10   0.07   0.14  0.01
## LanguageItalian                                       -0.20  -0.29  -0.12  0.00

2008-2010

m_GA_bam_2009 <- bam(GA_days ~ s(BW) +  s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) + parity_cat + sex + Urban3 + mother_nationality_cat2 + Language, data=bevn_eco_in7_2008_10)
summary(m_GA_bam_2009)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GA_days ~ s(BW) + s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) + 
##     parity_cat + sex + Urban3 + mother_nationality_cat2 + Language
## 
## Parametric coefficients:
##                                                       Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)                                          275.69818    0.03934 7008.938  < 2e-16 ***
## parity_cat(1,2]                                       -1.85794    0.03988  -46.593  < 2e-16 ***
## parity_cat(2,3]                                       -2.27191    0.06109  -37.189  < 2e-16 ***
## parity_cat(3,20]                                      -2.34282    0.10408  -22.510  < 2e-16 ***
## sex2                                                   1.81033    0.03583   50.528  < 2e-16 ***
## Urban3                                                 0.16255    0.03745    4.340 1.42e-05 ***
## mother_nationality_cat2Africa                         -0.17976    0.12735   -1.412  0.15810    
## mother_nationality_cat2Asia                           -0.87949    0.10418   -8.442  < 2e-16 ***
## mother_nationality_cat2Europe                         -0.29497    0.04100   -7.194 6.30e-13 ***
## mother_nationality_cat2Northern America               -0.58820    0.22410   -2.625  0.00867 ** 
## mother_nationality_cat2Oceania                        -0.14096    0.52659   -0.268  0.78894    
## mother_nationality_cat2Southern and Central America   -1.51662    0.12663  -11.976  < 2e-16 ***
## LanguageFrench                                        -0.08757    0.04398   -1.991  0.04649 *  
## LanguageItalian                                       -0.20814    0.09197   -2.263  0.02363 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df        F p-value    
## s(BW)            8.872  8.994 24566.37  <2e-16 ***
## s(mat_age)       7.177  7.939    66.93  <2e-16 ***
## s(birth_Y_M_num) 6.183  7.340    31.00  <2e-16 ***
## s(mean_ssep2)    4.779  5.919    20.22  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.505   Deviance explained = 50.5%
## fREML = 7.7699e+05  Scale est. = 68.838    n = 219790
  #plot(m_GA_bam_2009, all.terms = TRUE,pages=1, shade = TRUE, shade.col = "lightblue")
  plot(m_GA_bam_2009, select=1, shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_2009)[1]) #BW

  plot(m_GA_bam_2009, select=2, ylim=c(265, 280),shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_2009)[1]) #mat age

  plot(m_GA_bam_2009, select=3, ylim=c(274, 278),shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_2009)[1],  xlab  = "13 : 2008-01,   20: 2008-08,   30: 2009-06,   40: 2010-04,   50: 2011-02") #birth_Y_M_num

  plot(m_GA_bam_2009, select=4, ylim=c(272, 280),shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_2009)[1]) #ssep

  gam.check(m_GA_bam_2009)

## 
## Method: fREML   Optimizer: perf newton
## full convergence after 12 iterations.
## Gradient range [-0.000335078,0.0003341467]
## (score 776990.4 & scale 68.83792).
## Hessian positive definite, eigenvalue range [0.9452316,109886].
## Model rank =  50 / 50 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                    k'  edf k-index p-value   
## s(BW)            9.00 8.87    1.00   0.585   
## s(mat_age)       9.00 7.18    0.96   0.005 **
## s(birth_Y_M_num) 9.00 6.18    1.01   0.710   
## s(mean_ssep2)    9.00 4.78    1.00   0.475   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  concurvity(m_GA_bam_2009, full=TRUE)
##               para      s(BW) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst    0.7978125 0.04851831 0.13569927      0.005392540    0.18029328
## observed 0.7978125 0.01941720 0.11222580      0.004359862    0.07127053
## estimate 0.7978125 0.03114362 0.08700719      0.002902380    0.08241982
  #estimate, 95%CI and d
  beta <- coef(m_GA_bam_2009)
  Vb <- vcov(m_GA_bam_2009, unconditional = TRUE)
  n <- 219790
  se <- sqrt(diag(Vb))
  d <- (abs(beta)/(sqrt(n)*se))
  my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))
  round(head(my.ci,14), 2)
##                                                        beta    lci    uci     d
## (Intercept)                                          275.70 275.62 275.78 14.95
## parity_cat(1,2]                                       -1.86  -1.94  -1.78  0.10
## parity_cat(2,3]                                       -2.27  -2.39  -2.15  0.08
## parity_cat(3,20]                                      -2.34  -2.55  -2.14  0.05
## sex2                                                   1.81   1.74   1.88  0.11
## Urban3                                                 0.16   0.09   0.24  0.01
## mother_nationality_cat2Africa                         -0.18  -0.43   0.07  0.00
## mother_nationality_cat2Asia                           -0.88  -1.08  -0.68  0.02
## mother_nationality_cat2Europe                         -0.29  -0.38  -0.21  0.02
## mother_nationality_cat2Northern America               -0.59  -1.03  -0.15  0.01
## mother_nationality_cat2Oceania                        -0.14  -1.17   0.89  0.00
## mother_nationality_cat2Southern and Central America   -1.52  -1.76  -1.27  0.03
## LanguageFrench                                        -0.09  -0.17   0.00  0.00
## LanguageItalian                                       -0.21  -0.39  -0.03  0.00

2017-2020

m_GA_bam_covid <- bam(GA_days ~ s(BW) +  s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) + parity_cat + sex + Urban3 + mother_nationality_cat2 + Language, data=bevn_eco_in7_2017_20)
summary(m_GA_bam_covid)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## GA_days ~ s(BW) + s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) + 
##     parity_cat + sex + Urban3 + mother_nationality_cat2 + Language
## 
## Parametric coefficients:
##                                                       Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)                                          276.47736    0.03150 8777.856   <2e-16 ***
## parity_cat(1,2]                                       -2.04786    0.03143  -65.164   <2e-16 ***
## parity_cat(2,3]                                       -2.69664    0.04796  -56.228   <2e-16 ***
## parity_cat(3,20]                                      -2.85860    0.08356  -34.209   <2e-16 ***
## sex2                                                   1.80880    0.02837   63.757   <2e-16 ***
## Urban3                                                 0.05005    0.02973    1.683   0.0923 .  
## mother_nationality_cat2Africa                          0.91687    0.07848   11.682   <2e-16 ***
## mother_nationality_cat2Asia                           -0.94863    0.07386  -12.843   <2e-16 ***
## mother_nationality_cat2Europe                         -0.45937    0.03173  -14.477   <2e-16 ***
## mother_nationality_cat2Northern America               -0.38787    0.18813   -2.062   0.0392 *  
## mother_nationality_cat2Oceania                        -1.11572    0.46440   -2.402   0.0163 *  
## mother_nationality_cat2Southern and Central America   -1.84999    0.11179  -16.549   <2e-16 ***
## LanguageFrench                                         0.42080    0.03422   12.296   <2e-16 ***
## LanguageItalian                                       -0.11568    0.08240   -1.404   0.1604    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df         F p-value    
## s(BW)            8.956  8.999 38956.550 < 2e-16 ***
## s(mat_age)       7.036  7.817   164.828 < 2e-16 ***
## s(birth_Y_M_num) 3.077  3.826     4.686 0.00137 ** 
## s(mean_ssep2)    7.577  8.501    15.089 < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.517   Deviance explained = 51.7%
## fREML = 1.1677e+06  Scale est. = 65.374    n = 332752
  #plot(m_GA_bam_covid, all.terms = TRUE,pages=1, shade = TRUE, shade.col = "lightblue")
  plot(m_GA_bam_covid, select=1, shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_covid)[1]) #BW

  plot(m_GA_bam_covid, select=2, ylim=c(265, 280),shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_covid)[1]) #mat age

  plot(m_GA_bam_covid, select=3, ylim=c(274, 278),shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_covid)[1], xlab  = "120 : 2016-12,   130: 2017-10,   140: 2018-08,   150: 2019-06,   160: 2020-04") #birth_Y_M_num

  plot(m_GA_bam_covid, select=4, ylim=c(272, 280),shade = TRUE, shade.col = "lightblue",  shift = coef(m_GA_bam_covid)[1]) #ssep

  gam.check(m_GA_bam_covid)

## 
## Method: fREML   Optimizer: perf newton
## full convergence after 12 iterations.
## Gradient range [-0.0002086092,0.000205571]
## (score 1167708 & scale 65.37401).
## Hessian positive definite, eigenvalue range [0.6299583,166367].
## Model rank =  50 / 50 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                    k'  edf k-index p-value
## s(BW)            9.00 8.96    1.00    0.42
## s(mat_age)       9.00 7.04    1.01    0.82
## s(birth_Y_M_num) 9.00 3.08    1.00    0.50
## s(mean_ssep2)    9.00 7.58    1.01    0.70
  concurvity(m_GA_bam_covid, full=TRUE)
##               para      s(BW) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst    0.8021077 0.04825138  0.1070555     0.0012721330    0.17792321
## observed 0.8021077 0.01913284  0.0402550     0.0004357399    0.02758032
## estimate 0.8021077 0.02994665  0.0678886     0.0010076568    0.08028872
  #estimate, 95%CI and d
  beta <- coef(m_GA_bam_covid)
  Vb <- vcov(m_GA_bam_covid, unconditional = TRUE)
  n <- 332752
  se <- sqrt(diag(Vb))
  d <- (abs(beta)/(sqrt(n)*se))
  my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))
  round(head(my.ci,14), 2)
##                                                        beta    lci    uci     d
## (Intercept)                                          276.48 276.42 276.54 15.22
## parity_cat(1,2]                                       -2.05  -2.11  -1.99  0.11
## parity_cat(2,3]                                       -2.70  -2.79  -2.60  0.10
## parity_cat(3,20]                                      -2.86  -3.02  -2.69  0.06
## sex2                                                   1.81   1.75   1.86  0.11
## Urban3                                                 0.05  -0.01   0.11  0.00
## mother_nationality_cat2Africa                          0.92   0.76   1.07  0.02
## mother_nationality_cat2Asia                           -0.95  -1.09  -0.80  0.02
## mother_nationality_cat2Europe                         -0.46  -0.52  -0.40  0.03
## mother_nationality_cat2Northern America               -0.39  -0.76  -0.02  0.00
## mother_nationality_cat2Oceania                        -1.12  -2.03  -0.21  0.00
## mother_nationality_cat2Southern and Central America   -1.85  -2.07  -1.63  0.03
## LanguageFrench                                         0.42   0.35   0.49  0.02
## LanguageItalian                                       -0.12  -0.28   0.05  0.00

Pre-term birth: yes/no

2007-2020

m_PTB_bam<-bam(PTB~s(BW)+ s(mat_age) + s(birth_Y_M_num) + s(month) +s(mean_ssep2) + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7)
  summary(m_PTB_bam)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## PTB ~ s(BW) + s(mat_age) + s(birth_Y_M_num) + s(month) + s(mean_ssep2) + 
##     parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                      Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)                                          -4.27613    0.01550 -275.892  < 2e-16 ***
## parity_cat(1,2]                                       0.09992    0.01228    8.135 4.11e-16 ***
## parity_cat(2,3]                                       0.22338    0.01922   11.620  < 2e-16 ***
## parity_cat(3,20]                                      0.39947    0.03191   12.520  < 2e-16 ***
## Urban3                                               -0.05225    0.01142   -4.574 4.78e-06 ***
## LanguageFrench                                       -0.15425    0.01311  -11.770  < 2e-16 ***
## LanguageItalian                                      -0.19535    0.02835   -6.890 5.57e-12 ***
## mother_nationality_cat2Africa                        -0.13392    0.03577   -3.743 0.000182 ***
## mother_nationality_cat2Asia                          -0.10938    0.02933   -3.729 0.000192 ***
## mother_nationality_cat2Europe                         0.05210    0.01243    4.192 2.77e-05 ***
## mother_nationality_cat2Northern America               0.12821    0.07285    1.760 0.078403 .  
## mother_nationality_cat2Oceania                        0.25209    0.15753    1.600 0.109543    
## mother_nationality_cat2Southern and Central America   0.21445    0.03881    5.525 3.29e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df   Chi.sq p-value    
## s(BW)            8.345  8.863 97491.56 < 2e-16 ***
## s(mat_age)       6.424  7.331   146.40 < 2e-16 ***
## s(birth_Y_M_num) 3.072  3.822    80.01 < 2e-16 ***
## s(month)         5.348  6.475    15.94 0.01554 *  
## s(mean_ssep2)    5.483  6.674    23.79 0.00097 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.373   Deviance explained = 42.5%
## fREML = 1.5599e+06  Scale est. = 1         n = 1099328
    #adding intercept
#plot(m_PTB_bam, pages = 1, trans = plogis, shift = coef(m_PTB_bam)[1])
  plot(m_PTB_bam, trans = plogis, select=1,,shift = coef(m_PTB_bam)[1], shade = TRUE, shade.col = "lightblue") #BW

  plot(m_PTB_bam, trans = plogis, select=2, ylim=c(0, 0.05), shift = coef(m_PTB_bam)[1], shade = TRUE, shade.col = "lightblue") #matage

  plot(m_PTB_bam, trans = plogis, select=3, ylim=c(0.01, 0.02), shift = coef(m_PTB_bam)[1], shade = TRUE, shade.col = "lightblue",  xlab  = "1 : 2007-01,     25: 2009-01,      50: 2011-02,     75: 2013-03,     100: 2015-04,     125: 2017-05,     150: 2019-06") #birth_Y_M_num

  plot(m_PTB_bam, trans = plogis, select=4, ylim=c(0.012, 0.016), shift = coef(m_PTB_bam)[1], shade = TRUE, shade.col = "lightblue") #month

  plot(m_PTB_bam, trans = plogis, select=5, ylim=c(0.001, 0.02),shift = coef(m_PTB_bam)[1], shade = TRUE, shade.col = "lightblue") #ssep

  k.check(m_PTB_bam) #some kindex values are a bit far from 1
##                  k'      edf   k-index p-value
## s(BW)             9 8.345454 0.9139217  0.0525
## s(mat_age)        9 6.424215 0.9353811  0.4775
## s(birth_Y_M_num)  9 3.072089 0.9373408  0.5425
## s(month)          9 5.347513 0.9337870  0.4225
## s(mean_ssep2)     9 5.482747 0.9179691  0.1225
  concurvity(m_PTB_bam)
##             para      s(BW) s(mat_age) s(birth_Y_M_num)    s(month) s(mean_ssep2)
## worst    0.75695 0.02679970 0.12110615       0.03003623 0.028060770    0.17918840
## observed 0.75695 0.02174550 0.09399625       0.01560116 0.001701081    0.06553150
## estimate 0.75695 0.01814849 0.08072682       0.01101228 0.021035075    0.08951894
  #to calculate the OR of PTB
  OR <- exp(coef(m_PTB_bam))
  beta <- coef(m_PTB_bam)
  Vb <- vcov(m_PTB_bam, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 13), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.01 0.01 0.01 2.36
## parity_cat(1,2]                                      1.11 1.08 1.13 0.06
## parity_cat(2,3]                                      1.25 1.20 1.30 0.12
## parity_cat(3,20]                                     1.49 1.40 1.59 0.22
## Urban3                                               0.95 0.93 0.97 0.03
## LanguageFrench                                       0.86 0.84 0.88 0.09
## LanguageItalian                                      0.82 0.78 0.87 0.11
## mother_nationality_cat2Africa                        0.87 0.82 0.94 0.07
## mother_nationality_cat2Asia                          0.90 0.85 0.95 0.06
## mother_nationality_cat2Europe                        1.05 1.03 1.08 0.03
## mother_nationality_cat2Northern America              1.14 0.99 1.31 0.07
## mother_nationality_cat2Oceania                       1.29 0.94 1.75 0.14
## mother_nationality_cat2Southern and Central America  1.24 1.15 1.34 0.12

2008-2010

The variables ssep, month and birth_Y_M_num were set as non-smooth variables (model including them as smooth variables showed a linear relationship with PTB risk)

  m_PTB_bam_2009_2 <-bam(PTB~s(BW)+ s(mat_age) + birth_Y_M_num + month + mean_ssep2 + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2008_10)
  summary(m_PTB_bam_2009_2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## PTB ~ s(BW) + s(mat_age) + birth_Y_M_num + month + mean_ssep2 + 
##     parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                          -4.046320   0.098396 -41.123  < 2e-16 ***
## birth_Y_M_num                                        -0.005923   0.001215  -4.876 1.08e-06 ***
## month                                                 0.002883   0.003648   0.790 0.429447    
## mean_ssep2                                            0.001002   0.001454   0.689 0.490710    
## parity_cat(1,2]                                       0.041571   0.027001   1.540 0.123654    
## parity_cat(2,3]                                       0.180608   0.042228   4.277 1.89e-05 ***
## parity_cat(3,20]                                      0.338843   0.068907   4.917 8.77e-07 ***
## Urban3                                               -0.062921   0.024527  -2.565 0.010307 *  
## LanguageFrench                                       -0.124757   0.028516  -4.375 1.21e-05 ***
## LanguageItalian                                      -0.233616   0.058557  -3.990 6.62e-05 ***
## mother_nationality_cat2Africa                        -0.016162   0.088601  -0.182 0.855263    
## mother_nationality_cat2Asia                          -0.075412   0.068437  -1.102 0.270497    
## mother_nationality_cat2Europe                         0.081565   0.027526   2.963 0.003045 ** 
## mother_nationality_cat2Northern America               0.330363   0.148260   2.228 0.025863 *  
## mother_nationality_cat2Oceania                        0.481704   0.313962   1.534 0.124962    
## mother_nationality_cat2Southern and Central America   0.273003   0.081406   3.354 0.000798 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df   Chi.sq p-value    
## s(BW)      6.840  7.709 19924.46  <2e-16 ***
## s(mat_age) 5.576  6.569    42.04  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.362   Deviance explained = 41.2%
## fREML = 3.1193e+05  Scale est. = 1         n = 219790
  plot(m_PTB_bam_2009_2, trans = plogis, select=1, shade = TRUE, shade.col = "lightblue",  shift = coef(m_PTB_bam_2009_2)[1]) #BW

  plot(m_PTB_bam_2009_2, trans = plogis, select=2, ylim=c(0, 0.1),shade = TRUE, shade.col = "lightblue",  shift = coef(m_PTB_bam_2009_2)[1]) #mat age

  k.check(m_PTB_bam_2009_2)
##            k'      edf   k-index p-value
## s(BW)       9 6.839574 0.9350836  0.0325
## s(mat_age)  9 5.576486 0.9457077  0.1000
  concurvity(m_PTB_bam_2009_2)
##               para      s(BW) s(mat_age)
## worst    0.9850434 0.02716942 0.13455667
## observed 0.9850434 0.02211902 0.09906189
## estimate 0.9850434 0.01766117 0.08612003
  #to calculate the OR of PTB
  OR <- exp(coef(m_PTB_bam_2009_2))
  beta <- coef(m_PTB_bam_2009_2)
  Vb <- vcov(m_PTB_bam_2009_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 16), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.02 0.01 0.02 2.24
## birth_Y_M_num                                        0.99 0.99 1.00 0.00
## month                                                1.00 1.00 1.01 0.00
## mean_ssep2                                           1.00 1.00 1.00 0.00
## parity_cat(1,2]                                      1.04 0.99 1.10 0.02
## parity_cat(2,3]                                      1.20 1.10 1.30 0.10
## parity_cat(3,20]                                     1.40 1.23 1.61 0.19
## Urban3                                               0.94 0.89 0.99 0.03
## LanguageFrench                                       0.88 0.83 0.93 0.07
## LanguageItalian                                      0.79 0.71 0.89 0.13
## mother_nationality_cat2Africa                        0.98 0.83 1.17 0.01
## mother_nationality_cat2Asia                          0.93 0.81 1.06 0.04
## mother_nationality_cat2Europe                        1.08 1.03 1.15 0.05
## mother_nationality_cat2Northern America              1.39 1.04 1.86 0.18
## mother_nationality_cat2Oceania                       1.62 0.87 3.00 0.27
## mother_nationality_cat2Southern and Central America  1.31 1.12 1.54 0.15

2017-2020

#birth_Y_M, month and mean ssep as linear variables:
  m_PTB_bam_covid_2 <-bam(PTB~s(BW)+ s(mat_age) + birth_Y_M_num + month + mean_ssep2 + parity_cat + Urban3 + Language +  mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2017_20)
  summary(m_PTB_bam_covid_2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## PTB ~ s(BW) + s(mat_age) + birth_Y_M_num + month + mean_ssep2 + 
##     parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                          -4.3365644  0.1328509 -32.642  < 2e-16 ***
## birth_Y_M_num                                        -0.0009770  0.0007532  -1.297 0.194627    
## month                                                 0.0008431  0.0030337   0.278 0.781074    
## mean_ssep2                                            0.0006166  0.0012510   0.493 0.622082    
## parity_cat(1,2]                                       0.1638479  0.0228156   7.181  6.9e-13 ***
## parity_cat(2,3]                                       0.2983101  0.0355576   8.389  < 2e-16 ***
## parity_cat(3,20]                                      0.5173515  0.0576642   8.972  < 2e-16 ***
## Urban3                                               -0.0637814  0.0208789  -3.055 0.002252 ** 
## LanguageFrench                                       -0.2005760  0.0240309  -8.347  < 2e-16 ***
## LanguageItalian                                      -0.1866281  0.0561558  -3.323 0.000889 ***
## mother_nationality_cat2Africa                        -0.2119198  0.0622213  -3.406 0.000659 ***
## mother_nationality_cat2Asia                          -0.0952343  0.0519765  -1.832 0.066913 .  
## mother_nationality_cat2Europe                         0.0586032  0.0228280   2.567 0.010254 *  
## mother_nationality_cat2Northern America               0.0482582  0.1448537   0.333 0.739020    
## mother_nationality_cat2Oceania                        0.1284141  0.3413723   0.376 0.706790    
## mother_nationality_cat2Southern and Central America   0.2060507  0.0754333   2.732 0.006303 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df   Chi.sq  p-value    
## s(BW)      8.033  8.714 28326.69  < 2e-16 ***
## s(mat_age) 4.845  5.832    35.55 3.13e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.377   Deviance explained = 43.3%
## fREML = 4.7393e+05  Scale est. = 1         n = 332752
   plot(m_PTB_bam_covid_2, trans = plogis, select=1, shade = TRUE, shade.col = "lightblue",  shift = coef(m_PTB_bam_covid_2)[1]) #BW

  plot(m_PTB_bam_covid_2, trans = plogis, select=2, ylim=c(0, 0.05),shade = TRUE, shade.col = "lightblue",  shift = coef(m_PTB_bam_covid_2)[1]) #mat age

  k.check(m_PTB_bam_covid_2)
##            k'      edf   k-index p-value
## s(BW)       9 8.033448 0.9925878  0.9975
## s(mat_age)  9 4.844914 0.9445570  0.2325
  concurvity(m_PTB_bam_covid_2)
##               para      s(BW) s(mat_age)
## worst    0.9940987 0.02711559 0.10614541
## observed 0.9940987 0.02201811 0.04842693
## estimate 0.9940987 0.01778896 0.06719466
  #to calculate the OR of PTB
  OR <- exp(coef(m_PTB_bam_covid_2))
  beta <- coef(m_PTB_bam_covid_2)
  Vb <- vcov(m_PTB_bam_covid_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 16), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.01 0.01 0.02 2.40
## birth_Y_M_num                                        1.00 1.00 1.00 0.00
## month                                                1.00 0.99 1.01 0.00
## mean_ssep2                                           1.00 1.00 1.00 0.00
## parity_cat(1,2]                                      1.18 1.13 1.23 0.09
## parity_cat(2,3]                                      1.35 1.26 1.44 0.16
## parity_cat(3,20]                                     1.68 1.50 1.88 0.29
## Urban3                                               0.94 0.90 0.98 0.04
## LanguageFrench                                       0.82 0.78 0.86 0.11
## LanguageItalian                                      0.83 0.74 0.93 0.10
## mother_nationality_cat2Africa                        0.81 0.72 0.91 0.12
## mother_nationality_cat2Asia                          0.91 0.82 1.01 0.05
## mother_nationality_cat2Europe                        1.06 1.01 1.11 0.03
## mother_nationality_cat2Northern America              1.05 0.79 1.39 0.03
## mother_nationality_cat2Oceania                       1.14 0.58 2.22 0.07
## mother_nationality_cat2Southern and Central America  1.23 1.06 1.42 0.11

Low birth weight

2007-2020

Month variable was set as a linear term (when included as a smooth term, edf=1.3)

m_LBW_bam<-bam(LBW~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num) + month +s(mean_ssep2) + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7)
  summary(m_LBW_bam)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## LBW ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num) + month + s(mean_ssep2) + 
##     parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                        Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)                                          -4.1815402  0.0201764 -207.249  < 2e-16 ***
## month                                                 0.0006484  0.0017701    0.366   0.7141    
## parity_cat(1,2]                                      -0.6579649  0.0139870  -47.041  < 2e-16 ***
## parity_cat(2,3]                                      -0.7937344  0.0226928  -34.977  < 2e-16 ***
## parity_cat(3,20]                                     -0.8078470  0.0386216  -20.917  < 2e-16 ***
## Urban3                                                0.0172269  0.0125976    1.367   0.1715    
## LanguageFrench                                        0.2850660  0.0140620   20.272  < 2e-16 ***
## LanguageItalian                                       0.1779457  0.0302849    5.876 4.21e-09 ***
## mother_nationality_cat2Africa                         0.0007782  0.0376698    0.021   0.9835    
## mother_nationality_cat2Asia                          -0.0152786  0.0317249   -0.482   0.6301    
## mother_nationality_cat2Europe                        -0.1505381  0.0139476  -10.793  < 2e-16 ***
## mother_nationality_cat2Northern America              -0.4943755  0.0911825   -5.422 5.90e-08 ***
## mother_nationality_cat2Oceania                       -0.4153060  0.2025663   -2.050   0.0403 *  
## mother_nationality_cat2Southern and Central America  -0.3650949  0.0457429   -7.981 1.45e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df   Chi.sq p-value    
## s(GA_weeks)      7.933  8.460 89283.36  <2e-16 ***
## s(mat_age)       5.676  6.654   115.71  <2e-16 ***
## s(birth_Y_M_num) 3.911  4.836    55.59  <2e-16 ***
## s(mean_ssep2)    2.274  2.919   130.07  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.409   Deviance explained = 45.3%
## fREML = 1.5977e+06  Scale est. = 1         n = 1099328
    plot(m_LBW_bam, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_bam)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks

  plot(m_LBW_bam, trans = plogis, select=2, ylim=c(.005, .03), shift = coef(m_LBW_bam)[1], shade = TRUE, shade.col = "lightblue") #mat_age

  plot(m_LBW_bam, trans = plogis, select=3, ylim=c(.01, .02), shift = coef(m_LBW_bam)[1], shade = TRUE, shade.col = "lightblue",  xlab  = "1 : 2007-01,     25: 2009-01,      50: 2011-02,     75: 2013-03,    100: 2015-04,     125: 2017-05,    150: 2019-06") #birth_Y_M_num

  plot(m_LBW_bam, trans = plogis, select=4, ylim=c(.005, .025), shift = coef(m_LBW_bam)[1], shade = TRUE, shade.col = "lightblue") #ssep

  k.check(m_LBW_bam)
##                  k'      edf   k-index p-value
## s(GA_weeks)       9 7.933245 0.8882932  0.0000
## s(mat_age)        9 5.675613 0.9391889  0.1550
## s(birth_Y_M_num)  9 3.910952 0.9772231  0.9800
## s(mean_ssep2)     9 2.273894 0.9575589  0.5925
  concurvity(m_LBW_bam)
##               para s(GA_weeks) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst    0.8741049 0.016894030 0.12294860       0.03012542    0.17905289
## observed 0.8741049 0.007885243 0.04784953       0.01087336    0.12485367
## estimate 0.8741049 0.009156089 0.08218594       0.01143670    0.08949728
  #to calculate the OR of LBW for non smooth terms
  OR <- exp(coef(m_LBW_bam))
  beta <- coef(m_LBW_bam)
  Vb <- vcov(m_LBW_bam, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 14), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.02 0.01 0.02 2.31
## month                                                1.00 1.00 1.00 0.00
## parity_cat(1,2]                                      0.52 0.50 0.53 0.36
## parity_cat(2,3]                                      0.45 0.43 0.47 0.44
## parity_cat(3,20]                                     0.45 0.41 0.48 0.45
## Urban3                                               1.02 0.99 1.04 0.01
## LanguageFrench                                       1.33 1.29 1.37 0.16
## LanguageItalian                                      1.19 1.13 1.27 0.10
## mother_nationality_cat2Africa                        1.00 0.93 1.08 0.00
## mother_nationality_cat2Asia                          0.98 0.93 1.05 0.01
## mother_nationality_cat2Europe                        0.86 0.84 0.88 0.08
## mother_nationality_cat2Northern America              0.61 0.51 0.73 0.27
## mother_nationality_cat2Oceania                       0.66 0.44 0.98 0.23
## mother_nationality_cat2Southern and Central America  0.69 0.63 0.76 0.20

2008-2010

As the variables birth_Y_M_num and month are too correlated : I did the model without month variable. Also, birth_Y_M_num and mean_ssep2 were set as non smooth variables (edf ~1 fit hey were included as smooth variables)

m_LBW_bam_2009_2 <-bam(LBW~s(GA_weeks)+ s(mat_age) + birth_Y_M_num + mean_ssep2 + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2008_10)
  summary(m_LBW_bam_2009_2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## LBW ~ s(GA_weeks) + s(mat_age) + birth_Y_M_num + mean_ssep2 + 
##     parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                          -3.775443   0.109629 -34.438  < 2e-16 ***
## birth_Y_M_num                                         0.007496   0.001295   5.791 7.00e-09 ***
## mean_ssep2                                           -0.009640   0.001647  -5.853 4.82e-09 ***
## parity_cat(1,2]                                      -0.630716   0.031012 -20.338  < 2e-16 ***
## parity_cat(2,3]                                      -0.733919   0.050453 -14.547  < 2e-16 ***
## parity_cat(3,20]                                     -0.783168   0.085537  -9.156  < 2e-16 ***
## Urban3                                                0.052812   0.027693   1.907 0.056512 .  
## LanguageFrench                                        0.246728   0.031327   7.876 3.38e-15 ***
## LanguageItalian                                       0.296893   0.061089   4.860 1.17e-06 ***
## mother_nationality_cat2Africa                        -0.074333   0.095410  -0.779 0.435928    
## mother_nationality_cat2Asia                          -0.024941   0.075948  -0.328 0.742615    
## mother_nationality_cat2Europe                        -0.099987   0.031033  -3.222 0.001273 ** 
## mother_nationality_cat2Northern America              -0.725315   0.208002  -3.487 0.000488 ***
## mother_nationality_cat2Oceania                       -0.138900   0.388893  -0.357 0.720967    
## mother_nationality_cat2Southern and Central America  -0.423089   0.099577  -4.249 2.15e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df  Chi.sq p-value    
## s(GA_weeks) 7.110  7.829 18402.4 < 2e-16 ***
## s(mat_age)  3.831  4.745    21.7 0.00047 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.401   Deviance explained = 44.1%
## fREML = 3.1231e+05  Scale est. = 1         n = 219790
  plot(m_LBW_bam_2009_2, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks

  plot(m_LBW_bam_2009_2, trans = plogis, select=2, ylim=c(.005, .045), shift = coef(m_LBW_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #mat_age

  k.check(m_LBW_bam_2009_2)
##             k'      edf   k-index p-value
## s(GA_weeks)  9 7.109732 0.9524535  0.5325
## s(mat_age)   9 3.831235 0.9455602  0.3725
  concurvity(m_LBW_bam_2009_2)
##              para s(GA_weeks) s(mat_age)
## worst    0.984819 0.015475405 0.13702511
## observed 0.984819 0.007432975 0.03986371
## estimate 0.984819 0.008446321 0.08791625
  #to calculate the OR of LBW for non smooth terms
  OR <- exp(coef(m_LBW_bam_2009_2))
  beta <- coef(m_LBW_bam_2009_2)
  Vb <- vcov(m_LBW_bam_2009_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 15), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.02 0.02 0.03 2.09
## birth_Y_M_num                                        1.01 1.00 1.01 0.00
## mean_ssep2                                           0.99 0.99 0.99 0.01
## parity_cat(1,2]                                      0.53 0.50 0.57 0.35
## parity_cat(2,3]                                      0.48 0.43 0.53 0.41
## parity_cat(3,20]                                     0.46 0.39 0.54 0.43
## Urban3                                               1.05 1.00 1.11 0.03
## LanguageFrench                                       1.28 1.20 1.36 0.14
## LanguageItalian                                      1.35 1.19 1.52 0.16
## mother_nationality_cat2Africa                        0.93 0.77 1.12 0.04
## mother_nationality_cat2Asia                          0.98 0.84 1.13 0.01
## mother_nationality_cat2Europe                        0.90 0.85 0.96 0.06
## mother_nationality_cat2Northern America              0.48 0.32 0.73 0.40
## mother_nationality_cat2Oceania                       0.87 0.41 1.87 0.08
## mother_nationality_cat2Southern and Central America  0.66 0.54 0.80 0.23

2017-2020

Month variable was also not included in the model (correlation between birth_Y_M and month variables), and we set k=3 for s(birth_Y_M_num) because its almost linear (edf=1.99). Mean ssep set as linear term because almost linear (edf=1.47).

  m_LBW_bam_covid_2 <-bam(LBW~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num, k=3) + mean_ssep2 + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2017_20)
  summary(m_LBW_bam_covid_2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## LBW ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num, k = 3) + mean_ssep2 + 
##     parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                          -3.789475   0.086217 -43.953  < 2e-16 ***
## mean_ssep2                                           -0.008699   0.001378  -6.312 2.75e-10 ***
## parity_cat(1,2]                                      -0.692985   0.025756 -26.905  < 2e-16 ***
## parity_cat(2,3]                                      -0.864438   0.041814 -20.674  < 2e-16 ***
## parity_cat(3,20]                                     -0.804935   0.068636 -11.728  < 2e-16 ***
## Urban3                                                0.021859   0.023021   0.950  0.34236    
## LanguageFrench                                        0.293669   0.025656  11.446  < 2e-16 ***
## LanguageItalian                                       0.163459   0.059747   2.736  0.00622 ** 
## mother_nationality_cat2Africa                         0.072592   0.063694   1.140  0.25441    
## mother_nationality_cat2Asia                          -0.132195   0.057206  -2.311  0.02084 *  
## mother_nationality_cat2Europe                        -0.187461   0.025499  -7.352 1.96e-13 ***
## mother_nationality_cat2Northern America              -0.533030   0.172905  -3.083  0.00205 ** 
## mother_nationality_cat2Oceania                       -1.016028   0.490264  -2.072  0.03823 *  
## mother_nationality_cat2Southern and Central America  -0.263475   0.085099  -3.096  0.00196 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    Chi.sq p-value    
## s(GA_weeks)      7.293  7.922 26160.737  <2e-16 ***
## s(mat_age)       4.829  5.820    55.520  <2e-16 ***
## s(birth_Y_M_num) 1.705  1.913     2.865   0.163    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.41   Deviance explained = 46.1%
## fREML = 4.7418e+05  Scale est. = 1         n = 332752
  plot(m_LBW_bam_covid_2, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks

  plot(m_LBW_bam_covid_2, trans = plogis, select=2, ylim=c(.005, .06), shift = coef(m_LBW_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue") #mat_age

  plot(m_LBW_bam_covid_2, trans = plogis, select=3, ylim=c(.015, .025), shift = coef(m_LBW_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue", xlab  = "120 : 2016-12,   130: 2017-10,   140: 2018-08,   150: 2019-06,   160: 2020-04") #birth_Y_M_num

  k.check(m_LBW_bam_covid_2)
##                  k'      edf   k-index p-value
## s(GA_weeks)       9 7.292818 0.9667546  0.9375
## s(mat_age)        9 4.829053 0.9499728  0.5750
## s(birth_Y_M_num)  2 1.705205 0.9542470  0.7075
  concurvity(m_LBW_bam_covid_2)
##               para s(GA_weeks) s(mat_age) s(birth_Y_M_num)
## worst    0.9828781 0.018060995 0.10773748     0.0011409949
## observed 0.9828781 0.009161863 0.04775956     0.0007068171
## estimate 0.9828781 0.009997521 0.06895348     0.0010619952
  #to calculate the OR of LBW for non smooth terms
  OR <- exp(coef(m_LBW_bam_covid_2))
  beta <- coef(m_LBW_bam_covid_2)
  Vb <- vcov(m_LBW_bam_covid_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 14), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.02 0.02 0.03 2.09
## mean_ssep2                                           0.99 0.99 0.99 0.00
## parity_cat(1,2]                                      0.50 0.48 0.53 0.38
## parity_cat(2,3]                                      0.42 0.39 0.46 0.48
## parity_cat(3,20]                                     0.45 0.39 0.51 0.44
## Urban3                                               1.02 0.98 1.07 0.01
## LanguageFrench                                       1.34 1.28 1.41 0.16
## LanguageItalian                                      1.18 1.05 1.32 0.09
## mother_nationality_cat2Africa                        1.08 0.95 1.22 0.04
## mother_nationality_cat2Asia                          0.88 0.78 0.98 0.07
## mother_nationality_cat2Europe                        0.83 0.79 0.87 0.10
## mother_nationality_cat2Northern America              0.59 0.42 0.82 0.29
## mother_nationality_cat2Oceania                       0.36 0.14 0.95 0.56
## mother_nationality_cat2Southern and Central America  0.77 0.65 0.91 0.15

LBW vs “Normal” birthweight

2007-2020

Month was put as a linear variable.

m_LBW_normal_bam_2 <-bam(LBW_norm~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num) +s(mean_ssep2)+ month  + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7)
  summary(m_LBW_normal_bam_2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## LBW_norm ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) + 
##     month + parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                        Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)                                          -4.0369552  0.0196740 -205.192  < 2e-16 ***
## month                                                 0.0007164  0.0017715    0.404   0.6859    
## parity_cat(1,2]                                      -0.6409693  0.0140031  -45.773  < 2e-16 ***
## parity_cat(2,3]                                      -0.7685068  0.0227484  -33.783  < 2e-16 ***
## parity_cat(3,20]                                     -0.7778748  0.0387627  -20.068  < 2e-16 ***
## Urban3                                                0.0181157  0.0126123    1.436   0.1509    
## LanguageFrench                                        0.2804620  0.0140713   19.932  < 2e-16 ***
## LanguageItalian                                       0.1715774  0.0302897    5.665 1.47e-08 ***
## mother_nationality_cat2Africa                         0.0115914  0.0377714    0.307   0.7589    
## mother_nationality_cat2Asia                          -0.0143067  0.0317447   -0.451   0.6522    
## mother_nationality_cat2Europe                        -0.1445716  0.0139627  -10.354  < 2e-16 ***
## mother_nationality_cat2Northern America              -0.4899404  0.0912990   -5.366 8.04e-08 ***
## mother_nationality_cat2Oceania                       -0.4171044  0.2026645   -2.058   0.0396 *  
## mother_nationality_cat2Southern and Central America  -0.3583367  0.0458057   -7.823 5.16e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df   Chi.sq p-value    
## s(GA_weeks)      8.134  8.636 86595.52  <2e-16 ***
## s(mat_age)       5.636  6.616   117.38  <2e-16 ***
## s(birth_Y_M_num) 3.911  4.836    54.43  <2e-16 ***
## s(mean_ssep2)    2.325  2.984   132.31  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.408   Deviance explained = 44.5%
## fREML = 1.4478e+06  Scale est. = 1         n = 1011626
  plot(m_LBW_normal_bam_2, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_normal_bam_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks

  plot(m_LBW_normal_bam_2, trans = plogis, select=2, ylim=c(.005, .03), shift = coef(m_LBW_normal_bam_2)[1], shade = TRUE, shade.col = "lightblue") #mat_age

  plot(m_LBW_normal_bam_2, trans = plogis, select=3, ylim=c(.005, .025), shift = coef(m_LBW_normal_bam_2)[1], shade = TRUE, shade.col = "lightblue",  xlab  = "1 : 2007-01,     25: 2009-01,      50: 2011-02,     75: 2013-03,     100: 2015-04,    125: 2017-05,     150: 2019-06") #birth_Y_M_num

  plot(m_LBW_normal_bam_2, trans = plogis, select=4, ylim=c(.005, .030), shift = coef(m_LBW_normal_bam_2)[1], shade = TRUE, shade.col = "lightblue") #ssep

  k.check(m_LBW_normal_bam_2)
##                  k'      edf   k-index p-value
## s(GA_weeks)       9 8.133739 0.9273749  0.0225
## s(mat_age)        9 5.635822 0.9580069  0.4750
## s(birth_Y_M_num)  9 3.910756 0.9417167  0.1175
## s(mean_ssep2)     9 2.325108 0.9554838  0.3825
  concurvity(m_LBW_normal_bam_2)
##               para s(GA_weeks) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst    0.8728924 0.018337743 0.12132443       0.03011043    0.17987362
## observed 0.8728924 0.007720419 0.04750016       0.01120391    0.12435460
## estimate 0.8728924 0.009454512 0.07649425       0.01162304    0.08214456
  #to calculate the OR of LBW vs normal BW 
  OR <- exp(coef(m_LBW_normal_bam_2))
  beta <- coef(m_LBW_normal_bam_2)
  Vb <- vcov(m_LBW_normal_bam_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 14), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.02 0.02 0.02 2.23
## month                                                1.00 1.00 1.00 0.00
## parity_cat(1,2]                                      0.53 0.51 0.54 0.35
## parity_cat(2,3]                                      0.46 0.44 0.48 0.42
## parity_cat(3,20]                                     0.46 0.43 0.50 0.43
## Urban3                                               1.02 0.99 1.04 0.01
## LanguageFrench                                       1.32 1.29 1.36 0.15
## LanguageItalian                                      1.19 1.12 1.26 0.09
## mother_nationality_cat2Africa                        1.01 0.94 1.09 0.01
## mother_nationality_cat2Asia                          0.99 0.93 1.05 0.01
## mother_nationality_cat2Europe                        0.87 0.84 0.89 0.08
## mother_nationality_cat2Northern America              0.61 0.51 0.73 0.27
## mother_nationality_cat2Oceania                       0.66 0.44 0.98 0.23
## mother_nationality_cat2Southern and Central America  0.70 0.64 0.76 0.20

2008-2010

As always, month wasnt kept in the model (correlation between month and birth_Y_M_num variables). The variables birth_Y_M_num and ssep were set as linear terms.

m_LBW_normal_bam_2009_2 <-bam(LBW_norm~s(GA_weeks)+ s(mat_age) + birth_Y_M_num + mean_ssep2 + parity_cat +Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2008_10)
summary(m_LBW_normal_bam_2009_2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## LBW_norm ~ s(GA_weeks) + s(mat_age) + birth_Y_M_num + mean_ssep2 + 
##     parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                          -3.659351   0.109561 -33.400  < 2e-16 ***
## birth_Y_M_num                                         0.007519   0.001296   5.802 6.57e-09 ***
## mean_ssep2                                           -0.009721   0.001648  -5.898 3.68e-09 ***
## parity_cat(1,2]                                      -0.613393   0.031052 -19.754  < 2e-16 ***
## parity_cat(2,3]                                      -0.704519   0.050583 -13.928  < 2e-16 ***
## parity_cat(3,20]                                     -0.745120   0.085821  -8.682  < 2e-16 ***
## Urban3                                                0.055127   0.027710   1.989 0.046656 *  
## LanguageFrench                                        0.240118   0.031344   7.661 1.85e-14 ***
## LanguageItalian                                       0.294844   0.061125   4.824 1.41e-06 ***
## mother_nationality_cat2Africa                        -0.059372   0.095736  -0.620 0.535150    
## mother_nationality_cat2Asia                          -0.024569   0.076027  -0.323 0.746571    
## mother_nationality_cat2Europe                        -0.092698   0.031067  -2.984 0.002847 ** 
## mother_nationality_cat2Northern America              -0.725093   0.208405  -3.479 0.000503 ***
## mother_nationality_cat2Oceania                       -0.138224   0.388947  -0.355 0.722305    
## mother_nationality_cat2Southern and Central America  -0.414559   0.099779  -4.155 3.26e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df   Chi.sq  p-value    
## s(GA_weeks) 7.077  7.810 17825.65  < 2e-16 ***
## s(mat_age)  3.840  4.752    21.52 0.000506 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =    0.4   Deviance explained = 43.3%
## fREML = 2.8729e+05  Scale est. = 1         n = 201968
  plot(m_LBW_normal_bam_2009_2, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_normal_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks

  plot(m_LBW_normal_bam_2009_2, trans = plogis, select=2, ylim=c(.005, .05), shift = coef(m_LBW_normal_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #mat_age

  k.check(m_LBW_normal_bam_2009_2)
##             k'      edf   k-index p-value
## s(GA_weeks)  9 7.077019 0.9182415   0.010
## s(mat_age)   9 3.839971 0.9712118   0.965
  concurvity(m_LBW_normal_bam_2009_2)
##               para s(GA_weeks) s(mat_age)
## worst    0.9847922 0.016590747 0.13530354
## observed 0.9847922 0.006969109 0.04023792
## estimate 0.9847922 0.008342767 0.08762967
  #to calculate the OR of LBW vs normal BW 
  OR <- exp(coef(m_LBW_normal_bam_2009_2))
  beta <- coef(m_LBW_normal_bam_2009_2)
  Vb <- vcov(m_LBW_normal_bam_2009_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 15), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.03 0.02 0.03 2.02
## birth_Y_M_num                                        1.01 1.00 1.01 0.00
## mean_ssep2                                           0.99 0.99 0.99 0.01
## parity_cat(1,2]                                      0.54 0.51 0.58 0.34
## parity_cat(2,3]                                      0.49 0.45 0.55 0.39
## parity_cat(3,20]                                     0.47 0.40 0.56 0.41
## Urban3                                               1.06 1.00 1.12 0.03
## LanguageFrench                                       1.27 1.20 1.35 0.13
## LanguageItalian                                      1.34 1.19 1.51 0.16
## mother_nationality_cat2Africa                        0.94 0.78 1.14 0.03
## mother_nationality_cat2Asia                          0.98 0.84 1.13 0.01
## mother_nationality_cat2Europe                        0.91 0.86 0.97 0.05
## mother_nationality_cat2Northern America              0.48 0.32 0.73 0.40
## mother_nationality_cat2Oceania                       0.87 0.41 1.87 0.08
## mother_nationality_cat2Southern and Central America  0.66 0.54 0.80 0.23

2017-2020

Ssep variable was set as linear variables (when included as smooth variable, edf= 1.62)

  m_LBW_normal_bam_covid_2 <-bam(LBW_norm~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num) + mean_ssep2 + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2017_20)
summary(m_LBW_normal_bam_covid_2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## LBW_norm ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num) + mean_ssep2 + 
##     parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                          -3.618770   0.085835 -42.160  < 2e-16 ***
## mean_ssep2                                           -0.008824   0.001379  -6.400 1.55e-10 ***
## parity_cat(1,2]                                      -0.675156   0.025787 -26.182  < 2e-16 ***
## parity_cat(2,3]                                      -0.840120   0.041915 -20.044  < 2e-16 ***
## parity_cat(3,20]                                     -0.780183   0.068901 -11.323  < 2e-16 ***
## Urban3                                                0.022726   0.023047   0.986  0.32410    
## LanguageFrench                                        0.290190   0.025667  11.306  < 2e-16 ***
## LanguageItalian                                       0.153307   0.059736   2.566  0.01028 *  
## mother_nationality_cat2Africa                         0.081889   0.063870   1.282  0.19980    
## mother_nationality_cat2Asia                          -0.131873   0.057235  -2.304  0.02122 *  
## mother_nationality_cat2Europe                        -0.182846   0.025529  -7.162 7.94e-13 ***
## mother_nationality_cat2Northern America              -0.524457   0.173262  -3.027  0.00247 ** 
## mother_nationality_cat2Oceania                       -1.026005   0.489925  -2.094  0.03624 *  
## mother_nationality_cat2Southern and Central America  -0.256735   0.085221  -3.013  0.00259 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    Chi.sq p-value    
## s(GA_weeks)      7.269  7.907 25423.951  <2e-16 ***
## s(mat_age)       4.835  5.826    55.708  <2e-16 ***
## s(birth_Y_M_num) 1.856  2.316     3.334   0.262    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.409   Deviance explained = 45.4%
## fREML = 4.3553e+05  Scale est. = 1         n = 306277
 plot(m_LBW_normal_bam_covid_2, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_normal_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks

  plot(m_LBW_normal_bam_covid_2, trans = plogis, select=2, ylim=c(.005, .07), shift = coef(m_LBW_normal_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue") #mat_age

  plot(m_LBW_normal_bam_covid_2, trans = plogis, select=3, ylim=c(.020, .028), shift = coef(m_LBW_normal_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue", xlab  = "120 : 2016-12,   130: 2017-10,   140: 2018-08,   150: 2019-06,   160: 2020-04") #birth_Y_M_num

  k.check(m_LBW_normal_bam_covid_2)
##                  k'      edf   k-index p-value
## s(GA_weeks)       9 7.269236 0.9003570  0.0000
## s(mat_age)        9 4.835330 0.9640484  0.9425
## s(birth_Y_M_num)  9 1.855514 0.9436354  0.4875
  concurvity(m_LBW_normal_bam_covid_2)
##               para s(GA_weeks) s(mat_age) s(birth_Y_M_num)
## worst    0.9827989 0.020008547 0.10645213      0.001371803
## observed 0.9827989 0.009105054 0.04684491      0.001008478
## estimate 0.9827989 0.010278237 0.06782008      0.001042187
#to calculate the OR of LBW vs normal BW 
  OR <- exp(coef(m_LBW_normal_bam_covid_2))
  beta <- coef(m_LBW_normal_bam_covid_2)
  Vb <- vcov(m_LBW_normal_bam_covid_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 14), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.03 0.02 0.03 2.00
## mean_ssep2                                           0.99 0.99 0.99 0.00
## parity_cat(1,2]                                      0.51 0.48 0.54 0.37
## parity_cat(2,3]                                      0.43 0.40 0.47 0.46
## parity_cat(3,20]                                     0.46 0.40 0.52 0.43
## Urban3                                               1.02 0.98 1.07 0.01
## LanguageFrench                                       1.34 1.27 1.41 0.16
## LanguageItalian                                      1.17 1.04 1.31 0.08
## mother_nationality_cat2Africa                        1.09 0.96 1.23 0.05
## mother_nationality_cat2Asia                          0.88 0.78 0.98 0.07
## mother_nationality_cat2Europe                        0.83 0.79 0.88 0.10
## mother_nationality_cat2Northern America              0.59 0.42 0.83 0.29
## mother_nationality_cat2Oceania                       0.36 0.14 0.94 0.57
## mother_nationality_cat2Southern and Central America  0.77 0.65 0.91 0.14

Macrosomia vs “Normal” birthweight

2007-2020

m_macro_normal_bam <-bam(macro_norm~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num) + s(month) +s(mean_ssep2) + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7)
summary(m_macro_normal_bam)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## macro_norm ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num) + s(month) + 
##     s(mean_ssep2) + parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                       Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)                                          -2.988756   0.008479 -352.494  < 2e-16 ***
## parity_cat(1,2]                                       0.582219   0.008303   70.125  < 2e-16 ***
## parity_cat(2,3]                                       0.841384   0.011433   73.594  < 2e-16 ***
## parity_cat(3,20]                                      0.987128   0.018273   54.020  < 2e-16 ***
## Urban3                                                0.013847   0.007668    1.806  0.07095 .  
## LanguageFrench                                       -0.268400   0.009394  -28.571  < 2e-16 ***
## LanguageItalian                                      -0.372239   0.023024  -16.168  < 2e-16 ***
## mother_nationality_cat2Africa                         0.061706   0.021283    2.899  0.00374 ** 
## mother_nationality_cat2Asia                          -0.109246   0.022227   -4.915 8.88e-07 ***
## mother_nationality_cat2Europe                         0.171988   0.008118   21.186  < 2e-16 ***
## mother_nationality_cat2Northern America               0.337389   0.043411    7.772 7.73e-15 ***
## mother_nationality_cat2Oceania                        0.181090   0.104897    1.726  0.08428 .  
## mother_nationality_cat2Southern and Central America   0.148646   0.028587    5.200 1.99e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df    Chi.sq p-value    
## s(GA_weeks)      7.728  8.431 30433.750  <2e-16 ***
## s(mat_age)       4.525  5.449    12.595  0.0369 *  
## s(birth_Y_M_num) 4.301  5.295    77.959  <2e-16 ***
## s(month)         2.842  3.532     8.933  0.0468 *  
## s(mean_ssep2)    6.664  7.813    61.399  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0445   Deviance explained = 7.43%
## fREML = 1.4888e+06  Scale est. = 1         n = 1050380
plot(m_macro_normal_bam, trans = plogis, select=1, ylim=c(0, 0.6), shift = coef(m_macro_normal_bam)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks

  plot(m_macro_normal_bam, trans = plogis, select=2, ylim=c(0.03, 0.06), shift = coef(m_macro_normal_bam)[1], shade = TRUE, shade.col = "lightblue") #mat_age

  plot(m_macro_normal_bam, trans = plogis, select=3, ylim=c(0.045, 0.055), shift = coef(m_macro_normal_bam)[1], shade = TRUE, shade.col = "lightblue",  xlab  = "1 : 2007-01,     25: 2009-01,      50: 2011-02,     75: 2013-03,     100: 2015-04,    125: 2017-05,     150: 2019-06") #birth_Y_M_num

    plot(m_macro_normal_bam, trans = plogis, select=4, ylim=c(0.045, 0.052), shift = coef(m_macro_normal_bam)[1], shade = TRUE, shade.col = "lightblue") #month

  plot(m_macro_normal_bam, trans = plogis, select=5, ylim=c(0.02, 0.06), shift = coef(m_macro_normal_bam)[1], shade = TRUE, shade.col = "lightblue") #ssep

  k.check(m_macro_normal_bam)
##                  k'      edf   k-index p-value
## s(GA_weeks)       9 7.728394 0.9504278  0.8575
## s(mat_age)        9 4.525169 0.9498593  0.8175
## s(birth_Y_M_num)  9 4.301243 0.9400573  0.5875
## s(month)          9 2.842075 0.9521674  0.8775
## s(mean_ssep2)     9 6.664473 0.9509826  0.8350
  concurvity(m_macro_normal_bam)
##               para s(GA_weeks)  s(mat_age) s(birth_Y_M_num)     s(month) s(mean_ssep2)
## worst    0.7583176  0.01888534 0.125009429      0.030368293 0.0282971647    0.17894841
## observed 0.7583176  0.01351314 0.007078764      0.009518221 0.0008605803    0.02650765
## estimate 0.7583176  0.01216139 0.080052791      0.011622466 0.0211814067    0.08691995
#to calculate the OR of Macrosomia vs normal BW 
  OR <- exp(coef(m_macro_normal_bam))
  beta <- coef(m_macro_normal_bam)
  Vb <- vcov(m_macro_normal_bam, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 13), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.05 0.05 0.05 1.65
## parity_cat(1,2]                                      1.79 1.76 1.82 0.32
## parity_cat(2,3]                                      2.32 2.27 2.37 0.46
## parity_cat(3,20]                                     2.68 2.59 2.78 0.55
## Urban3                                               1.01 1.00 1.03 0.01
## LanguageFrench                                       0.76 0.75 0.78 0.15
## LanguageItalian                                      0.69 0.66 0.72 0.21
## mother_nationality_cat2Africa                        1.06 1.02 1.11 0.03
## mother_nationality_cat2Asia                          0.90 0.86 0.94 0.06
## mother_nationality_cat2Europe                        1.19 1.17 1.21 0.10
## mother_nationality_cat2Northern America              1.40 1.29 1.53 0.19
## mother_nationality_cat2Oceania                       1.20 0.98 1.47 0.10
## mother_nationality_cat2Southern and Central America  1.16 1.10 1.23 0.08

2008-2010

We set mat age and birth_Y_M_num variables as linear (edf= 1.005 and 1.378, respectively)

  m_macro_normal_bam_2009_2 <-bam(macro_norm~s(GA_weeks)+ s(mean_ssep2) + mat_age + birth_Y_M_num +parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2008_10)
  summary(m_macro_normal_bam_2009_2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## macro_norm ~ s(GA_weeks) + s(mean_ssep2) + mat_age + birth_Y_M_num + 
##     parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                                          -3.0624066  0.0594550 -51.508  < 2e-16 ***
## mat_age                                               0.0033146  0.0017206   1.926  0.05406 .  
## birth_Y_M_num                                        -0.0008767  0.0007783  -1.126  0.25997    
## parity_cat(1,2]                                       0.5519526  0.0184340  29.942  < 2e-16 ***
## parity_cat(2,3]                                       0.8393671  0.0252785  33.205  < 2e-16 ***
## parity_cat(3,20]                                      0.9752223  0.0400244  24.366  < 2e-16 ***
## Urban3                                                0.0148538  0.0168953   0.879  0.37931    
## LanguageFrench                                       -0.2539002  0.0210518 -12.061  < 2e-16 ***
## LanguageItalian                                      -0.3262404  0.0478543  -6.817 9.27e-12 ***
## mother_nationality_cat2Africa                         0.3452120  0.0536319   6.437 1.22e-10 ***
## mother_nationality_cat2Asia                          -0.1520186  0.0536794  -2.832  0.00463 ** 
## mother_nationality_cat2Europe                         0.1923115  0.0182076  10.562  < 2e-16 ***
## mother_nationality_cat2Northern America               0.2987413  0.0954085   3.131  0.00174 ** 
## mother_nationality_cat2Oceania                        0.1256268  0.2318870   0.542  0.58798    
## mother_nationality_cat2Southern and Central America   0.2611198  0.0577739   4.520 6.19e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                 edf Ref.df   Chi.sq p-value    
## s(GA_weeks)   6.602  7.499 6418.525  <2e-16 ***
## s(mean_ssep2) 2.339  3.010    6.209     0.1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.047   Deviance explained = 7.81%
## fREML = 2.9711e+05  Scale est. = 1         n = 209945
  plot(m_macro_normal_bam_2009_2, trans = plogis, select=1, ylim=c(0, 0.8), shift = coef(m_macro_normal_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks

    plot(m_macro_normal_bam_2009_2, trans = plogis, select=2, ylim=c(0.03, 0.06), shift = coef(m_macro_normal_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #ssep

  k.check(m_macro_normal_bam_2009_2)
##               k'     edf   k-index p-value
## s(GA_weeks)    9 6.60227 0.9345813  0.4100
## s(mean_ssep2)  9 2.33933 0.9479809  0.7875
  concurvity(m_macro_normal_bam_2009_2)
##               para s(GA_weeks) s(mean_ssep2)
## worst    0.9806332  0.01626321    0.17850048
## observed 0.9806332  0.01135410    0.16053032
## estimate 0.9806332  0.01127897    0.08817556
#to calculate the OR of Macrosomia vs normal BW 
  OR <- exp(coef(m_macro_normal_bam_2009_2))
  beta <- coef(m_macro_normal_bam_2009_2)
  Vb <- vcov(m_macro_normal_bam_2009_2, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 15), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.05 0.04 0.05 1.69
## mat_age                                              1.00 1.00 1.01 0.00
## birth_Y_M_num                                        1.00 1.00 1.00 0.00
## parity_cat(1,2]                                      1.74 1.68 1.80 0.30
## parity_cat(2,3]                                      2.31 2.20 2.43 0.46
## parity_cat(3,20]                                     2.65 2.45 2.87 0.54
## Urban3                                               1.01 0.98 1.05 0.01
## LanguageFrench                                       0.78 0.74 0.81 0.14
## LanguageItalian                                      0.72 0.66 0.79 0.18
## mother_nationality_cat2Africa                        1.41 1.27 1.57 0.19
## mother_nationality_cat2Asia                          0.86 0.77 0.95 0.08
## mother_nationality_cat2Europe                        1.21 1.17 1.26 0.11
## mother_nationality_cat2Northern America              1.35 1.12 1.63 0.17
## mother_nationality_cat2Oceania                       1.13 0.72 1.79 0.07
## mother_nationality_cat2Southern and Central America  1.30 1.16 1.45 0.14

2017-2020

m_macro_normal_bam_covid <-bam(macro_norm~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num) +s(mean_ssep2) + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2017_20)
summary(m_macro_normal_bam_covid)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## macro_norm ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) + 
##     parity_cat + Urban3 + Language + mother_nationality_cat2
## 
## Parametric coefficients:
##                                                       Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)                                          -2.966355   0.015301 -193.867  < 2e-16 ***
## parity_cat(1,2]                                       0.614600   0.015109   40.677  < 2e-16 ***
## parity_cat(2,3]                                       0.875458   0.020714   42.264  < 2e-16 ***
## parity_cat(3,20]                                      1.019343   0.033436   30.486  < 2e-16 ***
## Urban3                                                0.005643   0.013881    0.407 0.684368    
## LanguageFrench                                       -0.289867   0.016778  -17.277  < 2e-16 ***
## LanguageItalian                                      -0.421656   0.045138   -9.342  < 2e-16 ***
## mother_nationality_cat2Africa                        -0.051770   0.035678   -1.451 0.146767    
## mother_nationality_cat2Asia                          -0.105231   0.038140   -2.759 0.005796 ** 
## mother_nationality_cat2Europe                         0.130878   0.014712    8.896  < 2e-16 ***
## mother_nationality_cat2Northern America               0.299535   0.083080    3.605 0.000312 ***
## mother_nationality_cat2Oceania                       -0.058639   0.242115   -0.242 0.808630    
## mother_nationality_cat2Southern and Central America   0.116127   0.056445    2.057 0.039651 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df  Chi.sq  p-value    
## s(GA_weeks)      6.203  7.075 8864.53  < 2e-16 ***
## s(mat_age)       3.889  4.782   18.33  0.00175 ** 
## s(birth_Y_M_num) 2.444  3.044   12.81  0.00524 ** 
## s(mean_ssep2)    4.463  5.563   41.08 1.85e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.043   Deviance explained = 7.15%
## fREML = 4.5124e+05  Scale est. = 1         n = 318299
plot(m_macro_normal_bam_covid, trans = plogis, select=1, ylim=c(0, 0.7), shift = coef(m_macro_normal_bam_covid)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks

  plot(m_macro_normal_bam_covid, trans = plogis, select=2, ylim=c(0.02, 0.08), shift = coef(m_macro_normal_bam_covid)[1], shade = TRUE, shade.col = "lightblue") #mat_age

  plot(m_macro_normal_bam_covid, trans = plogis, select=3, ylim=c(0.04, 0.06), shift = coef(m_macro_normal_bam_covid)[1], shade = TRUE, shade.col = "lightblue", xlab  = "120 : 2016-12,   130: 2017-10,   140: 2018-08,   150: 2019-06,   160: 2020-04") #birth_Y_M_num

    plot(m_macro_normal_bam_covid, trans = plogis, select=4, ylim=c(0.02, 0.06), shift = coef(m_macro_normal_bam_covid)[1], shade = TRUE, shade.col = "lightblue") #ssep

  k.check(m_macro_normal_bam_covid)
##                  k'      edf   k-index p-value
## s(GA_weeks)       9 6.203012 0.9327167  0.1775
## s(mat_age)        9 3.888923 0.9544342  0.7525
## s(birth_Y_M_num)  9 2.444071 0.9590803  0.8325
## s(mean_ssep2)     9 4.463401 0.9251614  0.0475
  concurvity(m_macro_normal_bam_covid)
##               para s(GA_weeks) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst    0.7608652  0.02085818 0.11064626     0.0013996579    0.17772643
## observed 0.7608652  0.01556035 0.07319894     0.0008908752    0.01327785
## estimate 0.7608652  0.01460350 0.07032166     0.0009817982    0.07889043
#to calculate the OR of Macrosomia vs normal BW 
  OR <- exp(coef(m_macro_normal_bam_covid))
  beta <- coef(m_macro_normal_bam_covid)
  Vb <- vcov(m_macro_normal_bam_covid, unconditional = TRUE)
  se <- sqrt(diag(Vb))
  lci <- exp(beta-1.96*se)
  uci <- exp(beta+1.96*se)
  d <- abs(beta)/1.81
  my.ci <- data.frame(cbind(OR, lci, uci, d))
  round(head(my.ci, 13), 2)
##                                                        OR  lci  uci    d
## (Intercept)                                          0.05 0.05 0.05 1.64
## parity_cat(1,2]                                      1.85 1.79 1.90 0.34
## parity_cat(2,3]                                      2.40 2.30 2.50 0.48
## parity_cat(3,20]                                     2.77 2.60 2.96 0.56
## Urban3                                               1.01 0.98 1.03 0.00
## LanguageFrench                                       0.75 0.72 0.77 0.16
## LanguageItalian                                      0.66 0.60 0.72 0.23
## mother_nationality_cat2Africa                        0.95 0.89 1.02 0.03
## mother_nationality_cat2Asia                          0.90 0.84 0.97 0.06
## mother_nationality_cat2Europe                        1.14 1.11 1.17 0.07
## mother_nationality_cat2Northern America              1.35 1.15 1.59 0.17
## mother_nationality_cat2Oceania                       0.94 0.59 1.52 0.03
## mother_nationality_cat2Southern and Central America  1.12 1.01 1.25 0.06